Abstract
By utilizing spatiotemporal biological parameterizations, the adjoint variational method was applied to a 3D marine ecosystem dynamical model. The results of twin experiments demonstrated that the mean absolute error (MAE) of phytoplankton in the surface layer and the reduced cost function (RCF) could be used to evaluate both the simulation results and parameter estimation. Spatiotemporal variation of key parameters (KPs) was optimized in real experiments. The RCF and MAE in each assimilation period (72 periods per year) decreased obviously. The spatially varying KP (KPS), temporally varying KP (KPT), and constant KP (KPC) were obtained by averaging KPs of spatiotemporal variation. Another type of spatiotemporal KP (KPST) was represented by KPS, KPT, and KPC. The correlation analysis of KPs, either KPS or KPT, accorded with the real ecological mechanism. Running the model with KPS, KPT, KPC, and KPST, respectively, we found that MAE was the minimum when KPs were spatiotemporal variation (KPST), while MAE reached its maximum when KPs were constant (KPC). 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 a previous study). KPST, a representation of spatiotemporal variation, reduces the variable number in calculation.
1. 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–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–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 five-component 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–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 satellite-based 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 spatial parameterizations, which is a reference for estimating the spatially and temporally varying parameters in real experiments and future work. The marine ecosystem dynamical model and data including background field and observations are introduced in Section 2. Twin experiments are carried out, and the results are analyzed in Section 3. Based on the conclusion above, real experiments are implemented in Section 4. Conclusions are presented in the last section.
2. Model and Method
2.1. Ecosystem Model
The nitrogen-based ecological model built in this study contains four state variables: nutrient (), phytoplankton (), zooplankton (), and detritus (). 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 is the representation of the four state variables NPZD. The terms and represent the contribution to the concentration change of each state variable caused by biological processes and physical processes, respectively. The primitive equations of and are described in Appendix A. The parameters are listed in Table 1, and their initial values are mainly based on the previous studies [8, 15, 17, 26–28].
2.2. Data
2.2.1. Background Filed
The physical environment including background flow field and temperature field is described with monthly mean data from simple ocean data assimilation 1.4.2 (SODA). The idea behind SODA is to use direct observations to correct model errors in order to improve the reanalysis 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 horizontal grid;(3)the vertical levels are as follows: 5, 15, 25, 35, 46, 57, 70, 82, 96, 112, 129, 148, 171, 197, 229, 268, 317, 381, 465, 579, 729, 918, 1139, 1625, 1875, 2125, 2375, 2624, 2874, 3124, 3374, 3624, 3874, 4124, 4374, 4624, 4874, 5124, and 5374 m, the first fourteen of which are used for this study [29].
2.2.2. 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 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 and decreases rapidly when is greater than . 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 , is expressed as the form of (3) and calculated by using the SeaWiFS monthly mean data (averaged over 1998–2001) in January. 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 as an intermediate tool [31].
2.2.3. 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]: , where , and the half-saturation constant . 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: , where m, which is characteristic scale depth.
2.3. 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.
2.4. 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 . In this paper, is the assimilation step, is the direction (, 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:
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.
3. Twin Experiments and Analyses
Spatially varying parameters are estimated in this section. Based on the initial values of the ecological parameters listed in Table 1, two types of parameter variations are constructed: where 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.
3.1. 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 generated “observations” and are assimilated into the model to estimate the given spatial variation of .
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.
3.2. 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 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.
4. 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.
4.1. 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 the intermediate region, and other areas are denoted by C. KPs are constant in area C, and the values of KPs in area B are obtained from the values in area A and C by linear interpolation. So the KPs in area A are spatially varying 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 CF1 is about 3.9. After assimilation, all the CFs in the year decrease, and the maximum value of CF28 is less than 2.1. The RCF28, ratio of CF28 to CF1, 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 MAE1 is about 0.007, and the maximum value of MAE28 is less than 0.003. The ratio of MAE28 to MAE1 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.
(a)
(b)
(a)
(b)
4.2. 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.
(1) Spatially Varying KP. For each parameter of , we average the spatial variation fields in time domain and then obtain one field of spatial variation , . , , , , and are the spatial variations of , , , , and , which are shown in Figure 7.
(a)
(b)
(c)
(d)
(e)
(2) Temporally Varying KP. For each parameter of KP, we average parameter values in space, meaning that we obtain a constant KPT for each assimilation period and a time series for a whole year, , . KPT1, KPT2, KPT3, KPT4, KPT5 are the spatial variation of , , , , , which are shown in Figure 9.
(3) Constant KP. For each parameter of KP, average above time series in time domain (or average above spatial variation in space domain) and then we obtain a constant KPC, , . , KPC2, KPC3, KPC4, KPC5 are the constant , , , , and , the values of which are 0.5878, 0.4934, 0.0978, 0.054, and 0.0241, respectively.
(4) Spatially and Temporally Varying KP. For each parameter of KP, the spatial and temporal variations KPST are obtained by the following calculation: .
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 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 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.
5. 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.
Appendices
A. The Governing Equations of Biological Processes and Physical Processes in Biological Model
where , , and are the 3D water velocities; is the temperature of water; and are the horizontal and vertical eddy diffusivities, respectively. The parameters in the model are described in Table 1.
B. The Adjoint Model
where , , , and are the corresponding Lagrange adjoint operators of , , , and , respectively; is the available observations (phytoplankton) to be assimilated [17].
Acknowledgments
The authors would like to genuinely thank editor and reviewers for their hard work and the constructive suggestions to greatly improve this paper. Partial support for this research was provided by the National Natural Science Foundation of China through Grant nos. 41076006 and 41206001, Natural Science Foundation of Jiangsu Province of China through Grant no. BK2012315, the State Ministry of Science and Technology of China through Grant nos. 2013AA091201 and 2013AA091202, and the Fundamental Research Funds for the Central Universities 201261006 and 201262007.