The NorESM1 is a global fully coupled system, which is based on the Community Earth System Model version 1.0.3 (CESM1; Vertenstein et al., 2012). Unlike CESM1, NorESM1 uses the atmospheric component from the modified version of Community Atmosphere Model (CAM4-Oslo; Kirkevåg et al., 2013). The ocean physical component of NorESM1 is based on Miami Isopycnic Coordinate Ocean Model (MICOM; Bleck and Smith, 1990; Bleck et al., 1992) but with modified numerics and physics (Bentsen et al., 2012). The ocean BGC compartment in NorESM1 is the Hamburg Oceanic Carbon Cycle (HAMOCC; Maier-Reimer et al., 2005; Tjiputra et al., 2013), which is embedded with the isopycnic MICOM model (Assmann et al., 2010). The other components in the model are adopted in their original form from CESM1, which are the Community Land Model (CLM4; Oleson et al., 2010; Lawrence et al., 2011), the Los Alamos sea ice model (CICE4; Gent et al., 2011; Holland et al., 2012) and with the version 7 coupler (CPL7; Craig et al., 2012).

This study utilizes the medium-resolution version of NorESM1 (Tjiputra et al., 2013). The atmospheric component CAM4 and Land component CLM4 are configured on a horizontal resolution of 1.9° at latitude and 2.5° at longitude (approximately 2° finite volume grid). In the vertical, CAM4 features 26 hybrid sigma-pressure levels with model top at approximately 3 hPa. The ocean MICOM and sea ice CICE4 models have a common horizontal resolution of approximately 1° × 1° with refined grids near the Equator in meridional direction, and in both zonal and meridional direction at high latitudes. MICOM uses 51 isopycnal layers and 2 additional layers for representing the bulk mixed layer with time-evolving thicknesses and densities. The biogeochemical component HAMOCC utilizes the same spatial and temporal resolution as the ocean model.

HAMOCC includes an NPZD-type ecosystem module which was initially implemented by Six and Maier-Reimer (1996). It includes one generic class of phytoplankton, one generic class of zooplankton, three macronutrients (phosphate, nitrate, and silicate), and one micronutrient (dissolved iron). In addition to the ecosystem module, it also prognostically simulates full inorganic carbon chemistry, which includes dissolved inorganic carbon and alkalinity. Other key state variables include oxygen, dissolved organic carbon, particulate organic and inorganic carbon, and biogenic opal. The phytoplankton growth rate in the model is formulated as a function of temperature and light availability (Smith, 1936). The primary production is represented by a prognostic function of phytoplankton growth rate, which is limited by temperature, incoming shortwave radiation, and availability of nutrients. Further details of the HAMOCC can be sought from Tjiputra et al. (2013).

NorESM1 has been shown to capture the major observed modes of climatic variability (Bentsen et al., 2013). Further, many studies have demonstrated that it simulates well ENSO variability and its teleconnection (e.g., Sperber et al., 2013; Bellenger et al., 2014). In Anav et al. (2013), it was demonstrated that observed tropical inter-annual variability in ocean primary production reproduced by NorESM1 is in the top three among 18 ESM's used in their study. Tjiputra et al. (2013) evaluated the mean state of HAMOCC with NorESM and reported that NorESM satisfactorily reproduces many of the observed large scale ocean biogeochemical features.

The dual one step ahead smoother (DOSA, Gharamti et al., 2015) is an iterative smoother scheme based on the Ensemble Kalman Filter. The DOSA scheme respects the Bayesian consistency of the problem, and proceeds as shown in Figure 1. Here, we use the Deterministic EnKF (DEnKF) in DOSA. The DEnKF is a square-root (deterministic) formulation of the EnKF that solves the analysis without the need to perturb the observations. It inflates the errors by construction and is intended to perform well in operational applications (Sakov et al., 2012).

In the first step, the ensemble of analysed model state Xk−1a and its associated parameters at time k−1, Pk−1a, are integrated forward: Xkf=M(Xk-1a,Pk-1a). Note that the parameters are not changed during the model integration (i.e., Pk-1a=Pkf).

The observations at time k, y_k, are used to produce a smoothed estimate of the state and parameters at the previous analysis step k−1 as follows:

where,

where the superscript ^T denotes a matrix transpose and A the ensemble anomalies, i.e., A=X−x¯1T, with 1m=[1,1,…,1]∈ℛ1×N. In a similar way, the parameter ensemble, Pk−1s¯, at time k−1 is smoothed using y_k.

In the second step, the model is integrated forward to time k again but with smoothed ensemble of state Xk−1s and parameters Pk-1s; i.e., Xkf2=M(Xk-1s,Pk-1s).

The observations at time k, y^k, are then used again to produce an analysis of Xka

K_{k, k} is the standard Kalman gain and Akf2 are the ensemble anomalies constructed from Xkf2. It should be emphasized that the observations are used twice, but the second time they are used with a model that is using a different set of parameters. As such, model state is updated twice in a assimilation cycle but parameter is updated only at previous time step.

The model state (X) includes several ocean physical and biogeochemical prognostic model variables and they are updated in isopycnal coordinates as in Counillon et al. (2014); Wang et al. (2017). In the physical component, we update the full isopycnal temperature, salinity, layer thickness and velocities (53 isopycnal layers). Similarly, we update the biogeochemical variables at all isopycnal layers which include oxygen, phytoplankton, silicate, nitrate, total dissolved inorganic carbon, total alkalinity, dissolved organic carbon, particulate organic carbon, zooplankton, and biogenic silica. The list of selected BGC parameters for this study is provided in Section 2.3. For the assimilation system, the state vector is composed of the above physical and biogeochemical model variables along with the biogeochemical parameters. When one updates the ocean variable layer thickness, one effectively updates also the mass of the BGC quantities. In Bethke et al. (2021), it was shown that this approach conserves well BGC properties and does not introduce spurious upwelling at the Equator. However, with an EnKF, the linear analysis update returns unphysical values for non-Gaussian distributed variables. Some state variables have a physical constraint and their values should be positive definite such as layer thickness and tracer concentrations. For layer thickness we use the upscaling algorithm (Wang et al., 2016) while for the BGC concentration quantity (i.e., when updating the BGC state) a post-processing step is applied so that negative values are set to zero. We have not noticed degradations caused by the post-processing. We think that it is due to the fact that 1) part of the non-Gaussianity is already handled by the super layer algorithm which updates the layer thickness, and 2) with the smoothing flavor that the DOSA scheme provides, the updated parameters rarely became non-physical (i.e., negative).

Observations are used to update both ocean and BGC components in a strongly coupled framework (Penny and Hamill, 2017). The BGC component does not feedback to the physics in this version of NorESM and thus error in the physical state cannot be caused by error in the value of the BGC parameters. Therefore, for simultaneous state-parameter estimation, we update the parameter values from only BGC observations in Equations (1) and (3) while the state variables of ocean physics and BGC are updated using all available observations (see Figure 1).

The rest of the configuration for the assimilation experiments in this work follows that of Counillon et al. (2014, 2016); Wang et al. (2016, 2017) and Bethke et al. (2021). The assimilation algorithm uses a local analysis framework (Evensen, 2003; Sakov et al., 2012), where a local analysis is computed for one horizontal grid point at a time by utilizing all available observation in a spatial window around the grid point. A quasi-Gaussian and distance-dependent localization function (Gaspari and Cohn, 1999) is used to smooth the impact at the boundary of the localisation radius. In this work, the localization radius varies with latitude for both hydrographic profile and BGC observations (Wang et al., 2017). We do not use vertical localization. A moderation and a pre-screening technique (Sakov et al., 2012) is used to sustain the ensemble spread during the assimilation period. We also use the moderation technique, where observation error variance is increased (here by a factor of 4) for the update of the ensemble anomalies [Equation (3)] while the original value of the observation error variance is kept to update the ensemble mean [Equation (1)]. The pre-screening method inflates the observation error such that the analysis remains within two standard deviations of the forecast error from the ensemble mean.

We test the potential of the DOSA to optimise BGC parameters with NorCPM in an identical twin experiment framework. A reference model simulation performed with the prescribed parameter values is considered as the truth. We aim to retrieve the parameter values in the truth that are assumed to be unknown. We focus on optimizing five BGC parameters of NorESM1, which are among the most uncertain in the BGC model component. The parameters are: 1) the half-saturation constant for nutrient uptake during the phytoplankton growth (BKPHY), 2) Maximum zooplankton grazing rate (GRAZRA), 3) Phytoplankton exudation rate, i.e., the rate of dissolved organic carbon release by phytoplankton (GAMMAP), 4) Sinking speed for particulate organic carbon (WPOC), and 5) Half-saturation constant for silicate uptake during biogenic opal production (BKOPAL). A complete list of the ecosystem parameters used in the HAMOCC model is documented in Maier-Reimer et al. (2005).

The “true” parameters (TP) values are constant in time but they vary spatially (the first row of Figure 6). They have two Gaussian anomalies centered randomly: one with an isotropic distribution and another with an anisotropic distribution. The spatial pattern is purely artificial but more as a way to test the robustness of the proposed parameter estimation method in retrieving spatially varying pattern. The characteristic length scale and structure of those perturbations are unknown so that the DA system cannot be tuned specifically. The initial (first guess) perturbed parameters (PP) values (the second row of Figure 6) are sampled from a multi-variate Gaussian distribution, with a spatially uniform value for each ensemble member. The ensemble mean of PP is set intentionally to be 25% lower than the global mean of TP. The standard deviation of the ensemble of PP is equal to 33% of the global mean of TP. As such the PP are chosen so that the ensemble mean differs from the truth but that it encloses the truth value.

We have first performed a NorESM simulation (one member/realisation) with TP values (the first row of Figure 6) from 1980 to 1999 that henceforth referred to as TRUTH. It has been initialized in January 1980 from member one of the 30-member NorCPM1 historical simulation integrated with historical forcing from 1850 to 2014 following phase 6 of the Coupled Model Inter-comparison Project (CMIP6) protocol (Bethke et al., 2021). A tiny perturbation of 10^−6°C was added to SST of that member in January 1980. We constructed synthetic observations from monthly averages of TRUTH with random white noise taking into account for observational error. The observation error was specified to be equal to one standard deviation of the temporal variability in TRUTH. The observation error varies with grid cell and calendar month. The monthly averaged observations of temperature, salinity, phytoplankton concentration, Oxygen, Nitrate, Silicate, and Phosphate were chosen for this study. Synthetic observations of temperature and salinity were produced at 35 vertical z-levels sampling the full water column, while in horizontal direction we kept only points at every 5th model cell. The BGC observations have been produced at surface at every 5th grid cell. Observations in ice-covered water were discarded.

Four sets of experiments have been performed in this study as follows:

A summary of all experiments is given in Table 1.

A prerequisite for skillful parameter estimation is that assimilation constrains the error of the state variables well. Error in state variables, particularly at the surface and at inter-mediate depth can have a strong influence on near-surface biogeochemical processes. We work in a perfect model framework (i.e., all members use true parameter values) and we are interested in how well the error of the state variables is constrained by different observation networks. The monthly time evolution of the RMSE (averaged over the global domain and organized by depth level) of the NorESM free ensemble run (NorESM_TP) is shown in Figures 2, 3. We also show the RMSE-difference (RMSED) of EnKF_PHY, EnKF_BGC and EnKF_ALL assimilation experiments compared to NorESM_TP.

All assimilation experiments improve the accuracy of both the ocean and the BGC state variables in the near-surface levels compared to NorESM_TP. The difference between the three DA experiments is small. Sole assimilation of BGC or physical data alone is able to constrain well the surface ocean physical and biogeochemical variables and vice versa. This was somewhat unexpected, and it exemplifies well the potential of strongly coupled data assimilation (Penny and Hamill, 2017). It should be reminded that we assimilate in an isopycnal coordinate that has been shown to be more effective than assimilation in geopotential depth for surface observation (Gavart and De Mey, 1997; Counillon et al., 2016). The Analysis error for phytoplankton concentration is well reduced in all experiments. Similarly, clear improvements are shown for nutrients (phosphate and silicate) and oxygen estimates. BGC data assimilation alone yields the largest reduction of errors in the top 200 m. Below 200 m depth, the combined assimilation of physical and BGC observations provides slightly better performance than the other two experiments and it mitigates the degradation seen at deeper layers for some variables, e.g., temperature, phosphate, and oxygen. Overall, the accuracy of the combined assimilation experiment is slightly better. For example, the average error in EnKF_ALL for salinity in the top 1 km is 35% lower than that of the free run, while it is 24% for EnKF_PHY and 30% for EnKF_BGC. Similarly for Oxygen, EnKF_ALL has 28% lower error than the free run while EnKF_PHY is 19% and EnKF_BGC 24%.

This section presents and assesses spatially varying ecosystem parameters estimation with the DOSA scheme (reduction of error of the estimated parameter). We analyse the EnKF_PE, EnKF_PE_CO and EnKF_PE_SCO experiments (see Table 1). The time evolution of the global-averaged ensemble spreads and RMSEs of the estimated parameters obtained from EnKF_PE are shown in Figure 4. RMSEs for all parameters reduce with time and become stable within one year of assimilation. The reduction is largest for GRAZRA that shows a 72% error reduction from its initial distribution (Table 2). The error reduction in the remaining parameters WPOC, BKPHY, GAMMAP, and BKOPAL is about 54, 47, 40, and 30%, respectively. Similarly, the ensemble spreads of all parameters reduce with time. However, the reduction is quicker than for RMSE, which suggests that the system may benefit from using multiplicative or additive inflation (Mitchell and Houtekamer, 2000; Anderson, 2001). Similar results have been obtained from EnKF_PE_CO and EnKF_PE_SCO experiments (not shown).

Ensemble data assimilation estimates the parameters based on their correlation with the model misfits from observations (e.g., Anderson, 2001). In order to visualise the convergence process, we have used scatter plots (see Figure 5) of the parameter ensemble against the model deviation from the truth at different cycles of the assimilation experiment; i.e., in the start in January 1990 (red color), after the first assimilation, the second and the last assimilation cycles (green, purple, and blue colors, respectively). Each scattered dot represents one ensemble member and the big dots represent the ensemble mean. All ecosystem variables simulated using perturbed parameters show relatively large deviations from the TRUTH (y-axes in Figure 5), which underlines the sensitivity of the surface quantities to errors in the parameters. Among all parameters, GRAZRA shows the strongest linear relation with all variables and is the most important parameter for reducing model bias. For some parameters, the linear relationship is only strong with respect to some variables (e.g., WPOC with Silicate), which shows the importance of using multiple type of observations for the parameter estimation. The parameters are converging toward the true value very rapidly (already within the first 2 assimilation cycles), strongly reducing the error in ecosystem variables. After 11 assimilation cycles (corresponding to December 1990), one can notice that global means of the estimated parameters are very close to the true value and that the errors in the ecosystem variables are close to zero. This shows that the method converges quickly (within few assimilation steps) and is able to constrain the global mean of estimated parameters close to their true values by largely reducing the error in parameters. Similar results have been found for EnKF_PE_CO and EnKF_PE_SCO experiments (not shown).

We further analyse the spatial distribution of the estimated parameters obtained for December 1990 from all three experiments. Figure 6 shows the true and ensemble mean of the experiments with perturbed and estimated parameters. First, we note that the data assimilation yields a reduction of error compared to the initial values for all parameters. We can also notice that there is some spatial coherency in the value of the pattern retrieved. Parameters show relatively good agreement with the spatial distribution of the true value with spatial RMSE reduced by 75% for GRAZRA, by 50 % for BKPHY and WPOC, 40% for GAMMAP, and only 30% for BKOPAL, see Table 2. However, some differences remain and there is some small-scale noise. We suspect the latter to be related to spurious correlations present in our finite size ensemble (30). The places where the estimation fails to converge to true value may relate to places where the model is insensitive to the parameters. Hence, the parameter estimation can drift towards an erroneous value—e.g., as a response to spurious correlation or because of the approximation of Gaussianity and linearity during the analysis—without having an impact on the state error. In order to assess the impact of the error reduction on the state variable, we will freeze the parameter values at the last assimilation cycle and perform simulations in a free ensemble run and in a reanalysis mode in the following sections.

We verify the state accuracy of a free ensemble run that uses the final parameter's estimates of EnKF_PE, EnKF_PE_CO and EnKF_PE_SCO experiments, and refer to them as NorESM_PE, NorESM_PE_CO and NorESM_PE_SCO, respectively. We compare the performance against a free ensemble run using the true parameter (NorESM_TP) and one using the perturbed parameter (NorESM_PP).

In Figure 7, NorESM_PP poorly simulates the phytoplankton activity with a systematic overestimation of phytoplankton concentrations and generally longer spring blooms during the summer seasons compared to NorESM_TP. The runs that use the estimated parameters clearly outperform NorESM_PP, and perform closely to NorESM_TP. For instance, the global averaged RMSE of the phytoplankton concentration with NorESM_PE and NorESM_TP is roughly 51 and 54% lower than in NorESM_PP. Similarly, the biases in the parameter estimation experiments are significantly reduced to the level of that of TP. The performance of TP and PEs are also comparable in the upper 200-m as shown in Figure 7C. Below this depth, the model seems to have no sensitivity to these parameters and all experiments perform equivalently.

This conclusion is also verified for others BGC quantities such as silicate (Figure 8), phosphate (Figure 9), nitrate, and oxygen (not shown) in which PEs perform nearly as good as TP, while PP leads to a poor simulation of these variables and a large underestimate (negative bias) of the nutrient concentration, specifically during summer seasons, in the euphotic zone. This is related to the highly enhanced phytoplankton activity in PP, as seen earlier, which removes the available nutrients from the euphotic zone. RMSE and bias profile plots suggest that BGC parameters may strongly control the biogeochemical process in deeper layers up to roughly 500 m depth. Again, the simulated quality of nutrients and oxygen profiles by PEs and TP are very similar throughout entire depth.

We can also notice that NorESM_PE and NorESM_PE_SCO experiments show very comparable accuracy for ecosystem variables. Results from NorESM_PE_CO is not included to avoid overlapping lines in the Figures but we found very similar results. It suggests that BGC surface sparse climatological observations, and even very sparse climatological observations are somewhat sufficient to retrieve optimal ecosystem parameters with similar quality as time-varying observations using the DOSA-EnKF algorithm. Hence, in our model, the largest contribution to the error appears to be related to the seasonal cycle representation, which can be effectively corrected with a monthly climatology of observations.

We also analyse the impact of the parameters spatially. Figure 10 shows the RMSE and bias averaged over 100 m depth for the phytoplankton concentration and 500 m for the silicate profile. In general, bloom intensity increases from mid to high latitudes. It is seen that PP shows a larger overestimation of phytoplankton concentration, specifically over high latitudes and over some tropical regions, e.g., the eastern tropical Pacific Ocean. In all 3 PEs, the biases and RMSEs are reduced. The RMSE and bias patterns of the PEs closely match that of TP. In the case of silicate, PP shows a severe underestimation over the region where increased phytoplankton activities are seen.

Performance of the PEs is now assessed for an observation that was not used for tuning the parameters. Hence, we investigate air-sea CO₂ flux and net primary production (NPP) (Figure 11) and found consistent results. Similar to previous results, PEs maintain the accuracy closer to TP.

The above results provide strong evidence that the DOSA-EnKF system can successfully recover the optimal value for chosen BGC parameters. It also suggests that the remaining parametric error does not effectively influence the behavior of the model.

This section presents the accuracy of the state variables in the reanalyses which started in February 1990 and were run until December 1991 using fixed estimated parameters (REANA_PE). We compare the performance of REANA_PE with that of a reanalysis using perturbed parameters REANA_PP (lower benchmark) and a reanalysis using true parameters REANA_TP (upper benchmark).

The time evolution of the RMSE and bias over the euphotic zone is presented in Figure 12. First, we notice that assimilation improves the accuracy of the ecosystem variables in all experiments. For instance, the prior distributions of phytoplankton, oxygen and net primary production (NPP) are associated with high uncertainty and biases in the months with maximum bloom activity (January in the southern hemisphere). After a few assimilation cycles, there is a lower error during the bloom seasons. This can also be verified by comparing the phytoplankton free run accuracy shown in Figure 7 with the analysis accuracy presented here. For instance, the free run with PP shows RMSE values of roughly 70 μmolC m⁻³ for the spring bloom peak over the northern hemisphere (Figure 7), which is reduced to roughly 40 μmolC m⁻³ after assimilation (Figure 12).

Similar results are found for NPP (time evolution of free run not shown), for which observations are not being assimilated in the system (non-observed variable). Thus, the assimilation system is capable of improving the quality of not only the observed ecosystem variables (e.g., phytoplankton and Oxygen) but also of the non-observed variables. However, differences in the accuracy of reanalyses are clearly visible for the different set of parameters. REANA_PP shows larger uncertainty and stronger biases particularly during bloom seasons for phytoplankton, oxygen and net primary production reanalyses. The differences between REANA_PP and REANA_TP are more pronounced in the northern hemisphere than in the southern hemisphere. This exemplifies that assimilation cannot achieve optimum performance in the presence of model error.

The quality assessment of reanalysis has been further assessed in the deeper ocean by estimating the globally averaged RMSE and bias profiles for the top 1-km (Figure 13). The statistics have been computed using the July 1990 to December 1991 period. The first 6 months have been discarded for removing assimilation spinup that is longer in the deeper ocean. TP and PE show overall comparable performance for phytoplankton reanalysis whereas PP leads to degraded performance mostly in the top 150 m. Similar results can be seen for nutrients (e.g., phosphate and silicate) and oxygen profiles where differences of PP with TP or PE are more pronounced at inter-mediate depth levels. Still we can see that the performance of PE is not as efficient for oxygen and phosphate below 300 m. It would have been interesting to test whether training the parameters with deeper BGC observations (currently only available at the surface) would have improved performance there.