Optimization of the Spatiotemporal Parameters in a Dynamical Marine Ecosystem Model Based on the Adjoint Assimilation

1 Laboratory of Physical Oceanography, Ocean University of China, Qingdao 266100, China 2 Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China 3 Key Laboratory of Marine Spill Oil Identification and Damage Assessment Technology, The organaization of North China Sea Monitoring Center, Qingdao 266033, China 4Key Laboratory of Ocean Circulation and Waves, Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China


Introduction
Ecological dynamical models that are useful for filling in the gaps between satellite data can potentially provide more complete time and space distributions of the variables that observations cannot [1][2][3].But models are deficient in their representation of real ecological processes and interactions, and consequently their outputs stray from the observations.The combination of models and observations through data assimilation is an effective method that can provide a reasonable representation of biological variables to improve the estimation of a system's state.The errors associated with using either data or model alone are reduced in a complementary fashion [4].
Data assimilation with variational methods focuses on parameter estimation, which analyze and improve model performance by adjusting parameters to approach the observations [5].In contrast to modeling of physical ocean and atmosphere, the governing equations for biogeochemical models are complicated, and all the models contain a large number of parameters, many of which are poorly known.The model performance is highly dependent on those parameters whose values are difficult to be assigned, only a few of which are derived from observational data, and most of which are derived from laboratory studies [6].Data assimilation, in the parameter estimation context, is a systematic method in a way that is consistent with the available observations [7] and has been the focus of many researches over the past decade [8][9][10][11][12][13][14][15][16][17].
There are several techniques that can be used for parameter optimization (gradient descent methods, conjugate gradient algorithms, simulated annealing, microgenetic algorithm, Newton method, stochastic search algorithms), and the most common in ecosystem modeling is the adjoint method (essentially a gradient descent method) [3].It is an efficient optimization method that reduces the large number of required iterations by first computing the gradient of the cost function with respect to each parameter, although iteration is still required [3].Lawson et al. [18] were the first to assimilate data into a biogeochemical model by the variational adjoint method for recovering model parameters as well as initial conditions and then applied the same method to a fivecomponent model in a twin experiment.
Real marine ecosystem responds to changes in environmental conditions and biotic population, so the ecological parameters cannot strictly be regarded as constant in a study area especially on global scale.The parameter values tend to be specific to certain area and may not be applicable for others [8,[19][20][21], which is true for terrestrial ecosystems.There have already been several studies focusing on this problem.The model was optimized independently in Losa et al. 's [22] work, and adjusted parameters exhibited significant spatial variations.Losa et al. [22] then run a four-component marine ecosystem model using spatially varying parameters obtained from the above study with some smooth transition over domain boundaries and reproduced regional scale patterns in the SeaWiFS imagery.Hemmings et al. [11] also achieved an improvement in simulation by utilizing spatial variations in parameters.Doron et al. [23] used a 3D ocean coupled physical-biogeochemical model implemented on the North Atlantic at 1/4 ∘ and including six biogeochemical variables to estimate three parameters that were assumed to be stochastic and had regional variations by (1) Gaussian assumption and reduced rank approximation and (2) non-Gaussian and nonlinear extensions of the analysis step in twin experiments.Fan and Lv [17] assimilated SeaWiFS chlorophyll-a data into a simple NPZD model by the adjoint method in a climatological physical environment provided by FOAM.They selected five key parameters that were sensitive to the modeling status and uncorrelated with each other by sensitivity analysis and then explored a new method to recover parameters in which several grids were selected as independent grids, and the parameter values of other grids could be represented through linear interpolation of these independent grids.The feasibility of spatial parameterizations and the validity of adjoint model were justified.On the other hand, even if in the same study area, the physical environment also changes in response to seasonal variability in atmospheric conditions which determines the dynamics of the surface ocean mixed layer in which phytoplankton grows, and this can lead to changing proportion of the plankton species.Therefore, the parameters for ecological models should be time dependent.Mattern et al. [24] used a statistical emulator technique, the polynomial chaos expansion, to estimate time-dependent values for two parameters of a 3-dimensional biological ocean model.They obtained values for the phytoplankton carbon-to-chlorophyll ratio and the zooplankton grazing rate by minimizing the misfit between simulated and satellitebased surface chlorophyll.Solidoro et al. [25] presented the statistical analysis of seasonal and spatial variability of water quality parameters in the lagoon of Venice.
In this study, adjoint variational method is applied to an NPZD model.Our purpose is to find a relatively efficient and high-precision way of linear interpolation when adopting
Figure 1 shows a schematic diagram of the biological interactions among the four state variables.Dissolved nutrients are absorbed by phytoplankton, while phytoplankton is grazed by zooplankton.The mortality of phytoplankton and zooplankton is the main source of the detritus.The remineralization and sinking of the detritus are also considered in the biological processes.The primitive equations are presented in a united formula: where of ocean variables with a straightforward assimilation algorithm.The data information provided by SODA and used for this study such as velocity and temperature is introduced as follow: (1) the data spans the 44-year period from 1958 to 2001; (2) averages of model output variables are remapped onto a uniform global 0.5 ∘ × 0.5 ∘ horizontal grid; (3) the vertical levels are as follows:

Photosynthetically Active Radiation (PAR)
. Phytoplankton, the primary producer, is in charge of transforming inorganics into organics by photosynthesis in the upper layer in the ocean, and light is the essential condition for photosynthesis.The limitation factor for light in the equations is where  0 is the light intensity for maximum photosynthesis rate [30].This function represents that with the increasing of , photosynthesis rate increases exponentially when  is less than  0 and decreases rapidly when  is greater than  0 . par is the symbol of photosynthetically active radiation (PAR), a part of solar shortwave radiation.The wavelength of solar shortwave radiation absorbed by sea water is 100-1000 nm approximately, 49% of which is infrared radiation (longer than 780 nm) that is transformed into heat energy, 7% is ultraviolet radiation that is scattered by water in the upper ocean, and 44% is called PAR for photosynthesis.The solar shortwave radiation data used in this study is NCEP/NCAR reanalysis that uses a frozen state-of-the-art global data assimilation system and a database as complete as possible to achieve quality control and assimilate many sources of observations.The horizontal resolution is 1.875 ∘ × 1.875 ∘ , and hence, Cressman interpolation will be used in order to keep consistent with the resolution of the model in this study.We select a 10-year record averaged over the period of 1997-2006 at 6 hours intervals and calculate the climatological mean data.Diffuse attenuation coefficient for the photosynthetically available radiation in case-1 waters, denoted by   (PAR), is expressed as the form of (3) and calculated by using the SeaWiFS monthly mean data (averaged over 1998-2001) in January.  (490) is the diffuse attenuation coefficient for Case-1 waters when the wavelength of light is 490 nm and the adopted constant for pure water is 0.0166 m −1 , obtained by using [chl] as an intermediate tool [31].

Observation and Initial Field.
The phytoplankton data for model initialization and assimilation are converted from SeaWiFS monthly mean surface chlorophyll-a (chla) concentrations.To be consistent with the physical environment (climatological monthly mean data), the data which serve as climatological phytoplankton observations are averaged over the period of 1997-2002.
The state variable PNZD calculated in the model of this study is converted to nitrogen concentration, and the unit is mmol N m −3 .The nutrient data are WOA09 monthly mean nitrate data with a resolution of 1.0 ∘ latitude × 2.0 ∘ longitude.There are 15 vertical levels, the upper 10 of which are interpolated by linear interpolation in order to be in accordance with the model.When converting chla (mg m −3 ) to nitrate (mmol N m −3 ), we first convert the chla data to  (carbon) (mg m −3 ) through an empirical hyperbolic formula developed by Semovski and Wozniak [32]:  =  max (Chl  /(Chl  +  1/2 ))chl  , where  max = 90, and the half-saturation constant  1/2 = 0.477.Then a Redfield ratio (6.625) (Redfield et al. [33]) is adopted to convert  to nitrogen (mmol N m −3 ).
The biological model is initialized in January.The initial values of surface phytoplankton and the nitrate are mentioned above.Comparable climatologies do not exist for  and .Similar to Losa et al. [20], zooplankton and detritus in the surface layer are 0.02 and 0.1 mmol N m −3 , respectively.We assume that concentrations of , , and  decrease exponentially with depth in the initialization and take  as an example: () = (1) ⋅  −(()−(1))/ ch ( = 1, 2, . . ., 14), where  ch = 90 m, which is characteristic scale depth.

Model Setting.
The simulated domain in this study is 0.5 ∘ E-359.5 ∘ E, 74.5 ∘ S−88.5 ∘ N with a resolution of 1.0 ∘ latitude × 1.0 ∘ longitude.There are 14 vertical levels in the upper 200 m.The equations are solved by -grid in horizontal direction and -coordinate in vertical direction.All the data sets used in this model are interpolated to the simulation grids as described above.The time step is 6 hours.

Adjoint Method.
The adjoint model provides a method of calculating the gradient of cost function (CF) with respect to ecosystem parameter.A CF is usually used to measure the misfit.In this paper, the CF is defined in a weighted sum of squares: where  and  are the modeled surface phytoplankton values and the observations, respectively. refers to space while  refers to time.  is the spatial weight of the observation, which is defined as   =   /  , where   is the area of the grid, and   is the mean area of these grids which have observations.By introducing Lagrange multipliers, the adjoint model, equations of which are shown in Appendix B, can be constructed.The gradient is used to determine the direction to optimize the parameter because CF declines in the inverse direction of its gradient.Then the new values of the parameters, which are selected as control variables (KPs in this paper), can be obtained by KP +1 = KP  +     .In this paper,  is the assimilation step,   is the direction (  = (− ⋅ 5%/86400) ⋅ ⃗   /| ⃗   |, ⃗   is the gradient of cost function with respect to control variable), and   is step length, based on which the amount to modification of the control variables is made.The values of step-length are determined by inferring to a previous study [15], and the details are shown as follows, and the assimilation step is 28 for optimization: (5) 2.5.Spatial Biological Parameterizations.For details of spatial biological parameterizations, refer to previous study [17,34].At first, several grids are selected as independent grids, and the parameter values of these grids are independent.Further, an influencing radius  is selected to determine observations of how big an area (influencing zone) can be influenced by the parameter values of the independent grid. , is the parameter of independent grid (, ), and  , is that of grid (, ). , can be represented through linear interpolation of  , as the following relation [34]: where  ,,, is the coefficient of linear interpolation,  ,,, is the weighting factor, and  ,,, is the distance between the independent grid (, ) and the grid (, ).After  , is tuned,  , can be interpolated according to (6).The newly adjusted parameters are used to rerun the model.In this study, the distribution schemes of independent grids will be discussed in Section 3.1.Repeat the process for 28 times using the optimization scheme in Section 2.4.

Twin Experiments and Analyses
Spatially varying parameters are estimated in this section.
Based on the initial values CV 0 of the ecological parameters listed in Table 1, two types of parameter variations are constructed: Type 2: where CV is the parameter value scaled by the corresponding value listed in Table 1, and it is a function of latitude (74 ∘ S-90 ∘ N) and longitude (100 ∘ E-80 ∘ W).According to the two given types, the parameters under estimation show a variation of paraboloid of revolution type between 0.75 and 1.25 (scaled values).The experiment results are evaluated by comprehensive analysis of the MAE and RE of the estimated parameter.

Twin Experiment 1:
Distribution Schemes of the Independent Grids.Distribution schemes of the independent grids are discussed by estimating spatially varying   in this section.Give the spatial variation to   while other parameters in Table 1 remain unchanged as constants.The spatial variation is constructed as function (7).Run the model for 5 days with a 6-hour time step.The modeled data of phytoplankton in the surface layer are memorized as model The distribution schemes of independent grids are described as follows.
(A) There are 72 × 33 independent grids selected uniformly over the entire model domain, each from a 5 ∘ × 5 ∘ area, and  is 5 ∘ .
(B) There are 36 × 17 independent grids selected uniformly over the entire model domain, each from a 10 ∘ × 10 ∘ area, and  is 10 ∘ .
(C) There are 24 × 11 independent grids selected uniformly over the entire model domain, each from a 15 ∘ × 15 ∘ area, and  is 15 ∘ .
Results of the three twin experiments above are shown in Table 2.It is obvious that the more independent grids are used, the better the assimilation results are.

Twin Experiment 2: Estimating the Five Key Parameters (KPs).
The objective of this part is to verify the effectiveness of the model when it estimates more than one parameter at a time.  ,   , ,   , and   are more important than the other parameters and happen to constitute a set of uncorrelated parameters according to the sensitivity study in Fan and Lv' s work [17], so these five KPs are selected as control variables.Using the results from Section 3.1, we estimate the given spatial variations of KPs.Give the spatial variation type 1 to   ,   ,  and spatial variation type 2 to   ,   , run the model for five days, and the modeled data of phytoplankton in the surface layer are memorized as model generated "observations, " and estimate the five parameters simultaneously by assimilating these "observations." The MAE decreases to 0.0003 mmol N m −3 from 0.035 mmol N m −3 , and RCF decreases to 7.4 × 10 −4 .REs of KP are 0.37%, 0.50%, 1.09%, 0.65%, and 0.50%, respectively, most of which are less than 1%.The errors of KPs are bigger than 4% in Fan and Lv's [17] twin experiment, although the errors of the five parameters are under control.It indicates that the adjoint variational method in our model can present more satisfactory estimations of KPs.
The twin experiments above indicate that the adjoint variational method based on the spatial parameterizations is effective for estimating spatially varying parameters.In real experiments, we do not know the real distribution of ecological parameter before estimating, and usually only the cost function or MAE can be calculated to evaluate the assimilation performance.So we investigate the ability of MAE to assess simulation results.The relation between MAE and RE of parameters in twin experiments is shown by linear regression analysis in Figure 2. The solid line (regression line) is the result of linear regression in logarithmic coordinate system which can reflect the general trend of data.The function of regression line is  = 0.9 + 0.8 where a slope of 0.9, greater than zero, indicates that the relationship between RAE of parameters and MAE is direct ratio.Their correlation coefficient is 0.7, which indicates that MAE is reliable and can be used to evaluate the simulation and estimation.The correlation coefficient between MAE and RCF is 0.8 which indicates that MAE and RCF are useful in practical application because the value of ecological parameters is unknown and cannot be evaluation criteria for the model.

Real Experiments
The focus of most studies on parameter estimation is to analyze and improve model performance by adjusting parameters to conform to observations [8,26,35].In this section, our aim is to prove the rationality and necessity of KP's temporal and spatial variation by comparison in real experiments.

Real Experiment Based on Adjoint Assimilation.
Real experiment with adjoint variational method was carried out to optimize KP in this section by assimilating phytoplankton data in the surface layer which were converted from SeaWiFS ocean color data (see Section 2.2.3).The assimilation area (A) is between 17 ∘ N-45 ∘ N and 173 ∘ E-142 ∘ W as shown in Figure 3.A single assimilation period is 5-day long in this study (assimilating observational data for 5-day long).The assimilation processes are implemented one by one with the same first-guess values listed in Table 1.Our aim is to optimize KPs, with all the other parameters considered as constant during each assimilation period.KPs in area A are optimized based on the conclusion in Section 3.1, and the gray points in Figure 4 represent the independent grids.B is  for each assimilation period and are temporally varying for a whole year.This is also consistent with Garcia-Gorriz et al. 's [30] conclusion that different dynamics/periods produce different sets of optimal values of parameters.
Figure 5 shows the assimilation results.The minimum value of CF 1 is about 3.9.After assimilation, all the CFs in the year decrease, and the maximum value of CF 28 is less than 2.1.The RCF 28 , ratio of CF 28 to CF 1 , is smaller than 0.12 at most time of a year, meaning that the assimilation method is effective for all the cycles.The same conclusion can be obtained from the value of MAE.The minimum value of MAE 1 is about 0.007, and the maximum value of MAE 28 is less than 0.003.The ratio of MAE 28 to MAE 1 is smaller than 30 percent in a whole year.Figure 6 shows the results averaged over the 72 cycles.It is clear that the RCF descends rapidly at the first five assimilation steps, and then it reaches a steady state that is about 0.05.MAE also reaches a steady value of 0.002 after five assimilation steps, indicating that the variation of the parameters influences the model outputs a little.So we can say that the adjoint technique is computationally efficient to optimize the spatially and temporally varying parameters.

Comparative Experiments and Analysis.
After assimilation, we obtain the temporal and spatial variation of each KP, denoted by KP(, , ), where ,  refer to space while  refers to time.The values of KP on each grid point are time series for a whole year, and KP in each assimilation period is spatial variation field.In this section, spatially varying KP, temporally varying KP, constant KP, and spatially and temporally varying KP will be obtained, respectively, by the following way.
Obviously, the estimated KP shows great spatial and temporal variations in Figures 7 and 8, respectively.We find that   ,   , and  decrease while   and   increase and vice versa.So we divide the five parameters into two groups: (1)   ,   , and  and (2)   ,   .  is the maximum uptake rate of nutrient by phytoplankton.When   increases, the growth of phytoplankton will be faster.  is the mortality rate of zooplankton.When   increases, the amount of zooplankton will reduce; on the contrary, the amount of phytoplankton will increase. is the remineralization rate of detritus.When  increases, the remineralization will be accelerated to make more released nitrate available, which is propitious to the growth of phytoplankton.In a word, the bigger the three parameters are, the faster the phytoplankton grows.For the second group,   is zooplankton maximum )    1.
grazing rate, and   is the mortality rate of phytoplankton.When   decreases, the growth of phytoplankton becomes faster, and when   decreases, the amount of phytoplankton will increase.So the smaller the two parameters are, the faster the phytoplankton grows.The variations of KPs are consistent with each other in influencing phytoplankton, which accords with ecological mechanism in real ecosystem.
The same conclusion can be drawn from results of correlation analysis in Tables 3(a) and 3(b).The values in Table 3(a) are the correlation coefficients between any two spatial variations of KP, and the values in Table 3(b) are the correlation coefficients between any two temporal variations of KP.The correlation analysis shows that, in terms of either spatial variations or temporal variations, there was a positive correlation between any two of   ,   , and  in group (1) with the correlation coefficients reaching 0.93, and a positive correlation between   and   in group (2) with the correlation coefficients reaching 0.99, while there was a negative correlation between group (1) and group (2) with correlation coefficients reaching −0.99.
Run the model one year using KPs with KPS, KPT, KPC, and KPST, which are denoted by comparative experiment 1, comparative experiment 2, comparative experiment 3, and comparative experiment 4, respectively.The simulation results are compared with the real experiment in Section 4.1 obviously shown in Figure 9.After simulation based on adjoint method in real experiment, the MAEs in all periods are smaller than 0.004 mmol N m −3 , and the annual average value is about 0.0022 mmol N m −3 .By using KPST, the MAEs in all periods are also smaller than 0.004 mmol N m −3 , and the annual MAE is about 0.0029 mmol N m −3 , which is close to the results above.The annual average value of MAE is 0.007 mmol N m −3 when considering only temporally varying KPT and 0.008 mmol N m −3 when considering only spatially varying KPS, which is the scheme in Fan and Lv's [17] work.The MAE is larger than 0.01 mmol N m −3 , and the annual average value of MAE is 0.04 mmol N m −3 , which is 20 times larger than 0.002 mmol N m −3 (the results of real experiment and comparative experiment 4) when KPs are considered as a set of constants to run the model, which is the form of parameters in previous study [26,35].So we draw the conclusion that in marine ecosystem dynamical model, using KPs with spatiotemporal variation is reasonable and consistent with real ecosystems.Using spatiotemporal KPs could improve simulation precision compared with only using spatially varying KPs, temporally varying KPs, or constant KPs (these forms are the results in previous study).On the other hand, an arbitrary KP(, , ) with spatiotemporal variation, which is a three-dimensional array, can be represented by a product of two-dimensional array and one-dimensional array, so it will reduce the variable number for long-time model calculation.

Conclusion
The adjoint variational method is applied to an NPZD marine ecosystem dynamics model in this study.In twin experiments the modeled data of phytoplankton in the surface layer are memorized as model generated "observations" that are assimilated into the model to estimate the given spatial variations of parameters.We obtain the optimal scheme of independent points and find that the more independent grids we use, the better the assimilation results are.When estimating 5 KPs simultaneously, the given spatial variations of parameters can be successfully inverted with the RE of parameters less than 2%, by which the effectiveness of adjoint variational method based on spatial parameterizations is verified.All the results in twin experiments are analyzed by linear regression method, and we can draw the conclusion that the relationship between RE of parameters and MAE is of direct ratio with a correlation coefficient of 0.7, and the relationship between RCF and MAE is also of direct ratio with a correlation coefficient of 0.8.So the MAE and RCF can be used to evaluate the simulation and estimation, which is useful in practical application.
In real experiments, the feasibility of spatial and temporal parameterizations is tested.We take a square area in north pacific as a study area, and the climatological monthly mean SeaWiFS data in this area are assimilated into the model, with each assimilation period being 5 days and 72 periods in a whole year.We draw the conclusion that adjoint assimilation is a useful tool to optimize spatial and temporal variation of KPs.The RCF averaged over the 72 periods reduces to 5%, and MAE averaged over the 72 periods reduces from 0.015 mmol N m −3 to 0.002 mmol N m −3 after assimilation.
The correlation analysis of KPs, either spatially or temporally, accords with the ecological mechanism in real ecosystem, so the 5 KPs can be classified into two groups by the correlation.We obtain the KPS, KPT, KPC, and KPST, respectively.Run the model by using KPs with the above 4 kinds of variation, and compare the results with real experiment.We find that MAE is the minimum when KPs are spatiotemporal variation (KPST), and MAE is the maximum when KPs are constant.Therefore, considering parameters as constants, only spatial variation or only temporal variation (these forms are the results in a previous study) in parameter estimation and marine ecosystem modeling on large space and time scales is doubtable.In addition, it reduces the variable number in model calculation when a three-dimensional array is represented by the product of a two-dimensional array and a one-dimensional array.
It is demonstrated that the adjoint assimilation method based on spatial and temporal parameterizations is an effective tool to estimate parameters in the marine ecosystem.The simulation precision and efficiency are promoted by utilizing more independent grids.The improvements achieved in this study are helpful for research of numerical simulation in marine ecosystem dynamics in the future.The spatiotemporal variation of parameters, which is the focus of several researchers, is reasonable and necessary.It provides reference for future study and needs further research.The model improved by using spatiotemporal parameters can be coupled with global ocean-atmosphere model to study the global climate change.The method of estimating spatiotemporal parameters can be used in other numerical simulation.where  * ,  * ,  * , and  * are the corresponding Lagrange adjoint operators of , , , and , respectively;  is the available observations (phytoplankton) to be assimilated [17].

Figure 2 :
Figure 2: Statistical analysis in twin experiments.The solid line is regression line, and the space between short dash line is the 95% confidence interval.

Figure 3 :
Figure 3: Scope of assimilation area (noted by A).

Figure 5 :Figure 6 :
Figure 5: Assimilation results in Section 4.1.The horizontal coordinate indicates the days; the vertical coordinate represents CF on the left hand and MAE on the right hand.Before assimilation, the values of CF and MAE are denoted by CF 1 and MAE 1 ; after assimilation, the values of CF and MAE are denoted by CF 28 and MAE 28 .The statistical area is A as shown in Figure 3.

Figure 7 :
Figure 7: Spatial variation of KPs averaged in time domain.The values are scaled by the corresponding parameter values listed in Table1.

− 3 )Figure 9 :
Figure9: Comparison of MAE in logarithmic coordinates.Red line denotes result using only spatially varying KPS; blue dotted line denotes result using only temporally varying KPT; black dotted line denotes result using constant KPC; blue solid line denotes result using KP based on adjoint assimilation; and black solid line denotes result using KPST.

Table 1 :
is the representation of the four state variables NPZD.The terms Bio  and Phy  represent the contribution to the concentration change of each state variable caused by biological processes and physical processes, respectively.Ecological parameters and their initial values in the model.

Table 2 :
Comparison of the three strategies.generated "observations" and are assimilated into the model to estimate the given spatial variation of   .

Table 3 :
Time series of KPs averaged in space.The values are scaled by the corresponding parameter values listed in Table 1.(a) Correlation coefficients between any two parameters of spatially varying KP.(b) Correlation coefficients between any two parameters of temporally varying KP.