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

Authors: Xiaoyan Li; Chunhui Wang; Wei Fan; Xianqing Lv
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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 Section2. 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 (N), phytoplankton (P), zooplankton (Z), and detritus (D). 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:
(1)∂Ci∂t=Phyi+Bioi,(i=P,N,Z,D),
where Ci is the representation of the four state variables NPZD. The terms Bioi and Phyi represent the contribution to the concentration change of each state variable caused by biological processes and physical processes, respectively. The primitive equations of Bioi and Phyi 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].Table 1
Ecological parameters and their initial values in the model.
Symbol
Parameter
Value
Unit
V
m
Maximum uptake rate of nutrient by phytoplankton
0.8
Day−1
K
s
Half-saturation constant for nutrient uptake
0.25
(N) mmol· m−3
A
Q
10
Temperature coefficient for phytoplankton growth at 10°C
2.08
—
W
P
Sinking velocity of phytoplankton
1.0
m · day−1
G
m
Zooplankton maximum grazing rate
0.50
day−1
λ
Zooplankton Ivlev constant
0.20
(N)
B
Q
10
Temperature coefficient for zooplankton growthat 10°C
3.10
—
D
z
Zooplankton mortality rate
0.05
Day−1
γ
Assimilation ratio for zooplankton
0.30
—
D
p
Phytoplankton mortality rate
0.10
Day−1
e
Remineralization rate of detritus
0.0212
Day−1
θ
Excretion ratio for zooplankton
0.40
—
W
d
Sinking velocity of detritus
1.00
m · day−1Figure 1
Diagram of the compartments in the ecosystem model.


### 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 global0.5°×0.5° 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(2)f(I,z)=II0·exp(1-II0),I=Ipar·exp(Kpar·z),
where I0 is the light intensity for maximum photosynthesis rate [30]. This function represents that with the increasing of I, photosynthesis rate increases exponentially when I is less than I0 and decreases rapidly when I is greater than I0. Ipar 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 byKd(PAR), is expressed as the form of (3) and calculated by using the SeaWiFS monthly mean data (averaged over 1998–2001) in January. Kd(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].
(3)Kd(PAR)=0.0864+0.884Kd(490)-0.00137[Kd(490)]-1,Kd(490)=0.0166+0.0773[chl]0.6715.


#### 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 C (carbon) (mg m−3) through an empirical hyperbolic formula developed by Semovski and Wozniak [32]: C=ρmax(Chla/(Chla+K1/2))chla, where ρmax=90, and the half-saturation constant K1/2=0.477. Then a Redfield ratio (6.625) (Redfield et al. [33]) is adopted to convert C 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 forZ and D. 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 P, Z, and D decrease exponentially with depth in the initialization and take P as an example: P(k)=P(1)·e-(z(k)-z(1))/zch(k=1,2,…,14), where zch=90 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 byA-grid in horizontal direction and Z-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:(4)CF=12∑s,tWs(xs,t-x-s,t)2,
where x and x- are the modeled surface phytoplankton values and the observations, respectively. s refers to space while t refers to time. Ws is the spatial weight of the observation, which is defined as Ws=As/A-s, where As is the area of the grid, and A-s is the mean area of these grids which have observations.By introducing Lagrange multipliers, the adjoint model, equations of which are shown in AppendixB, 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 KPk+1=KPk+αkdk. In this paper, k is the assimilation step, dk is the direction (dk=(-V·5%/86400)·Gk→/|Gk→|, G→k is the gradient of cost function with respect to control variable), and αk 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)αk=1-0.1·(k-1),k=1,2,3,…,10,αk=0.1-0.01·(k-10),k=11,12,…,19,αk=0.01-0.001·(k-19),k=20,21,…,28.


### 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 R is selected to determine observations of how big an area (influencing zone) can be influenced by the parameter values of the independent grid. Eii,jj is the parameter of independent grid (ii, jj), and Fi,j is that of grid (i, j). Fi,j can be represented through linear interpolation of Eii,jj as the following relation [34]:
(6)Fi,j=∑ii,jjϕi,j,ii,jj·Eii,jj,ϕi,j,ii,jj=wi,j,ii,jj∑ii,jjwi,j,ii,jj,wi,j,ii,jj=R2-ri,j,ii,jj2R2+ri,j,ii,jj2,
where ϕi,j,ii,jj is the coefficient of linear interpolation, wi,j,ii,jj is the weighting factor, and ri,j,ii,jj is the distance between the independent grid (ii, jj) and the grid (i, j). After Eii,jj is tuned, Fi,j 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.


## 2.1. Ecosystem Model

The nitrogen-based ecological model built in this study contains four state variables: nutrient (N), phytoplankton (P), zooplankton (Z), and detritus (D). 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:
(1)∂Ci∂t=Phyi+Bioi,(i=P,N,Z,D),
where Ci is the representation of the four state variables NPZD. The terms Bioi and Phyi represent the contribution to the concentration change of each state variable caused by biological processes and physical processes, respectively. The primitive equations of Bioi and Phyi 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].Table 1
Ecological parameters and their initial values in the model.
Symbol
Parameter
Value
Unit
V
m
Maximum uptake rate of nutrient by phytoplankton
0.8
Day−1
K
s
Half-saturation constant for nutrient uptake
0.25
(N) mmol· m−3
A
Q
10
Temperature coefficient for phytoplankton growth at 10°C
2.08
—
W
P
Sinking velocity of phytoplankton
1.0
m · day−1
G
m
Zooplankton maximum grazing rate
0.50
day−1
λ
Zooplankton Ivlev constant
0.20
(N)
B
Q
10
Temperature coefficient for zooplankton growthat 10°C
3.10
—
D
z
Zooplankton mortality rate
0.05
Day−1
γ
Assimilation ratio for zooplankton
0.30
—
D
p
Phytoplankton mortality rate
0.10
Day−1
e
Remineralization rate of detritus
0.0212
Day−1
θ
Excretion ratio for zooplankton
0.40
—
W
d
Sinking velocity of detritus
1.00
m · day−1Figure 1
Diagram of the compartments in the ecosystem model.


## 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 global0.5°×0.5° 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(2)f(I,z)=II0·exp(1-II0),I=Ipar·exp(Kpar·z),
where I0 is the light intensity for maximum photosynthesis rate [30]. This function represents that with the increasing of I, photosynthesis rate increases exponentially when I is less than I0 and decreases rapidly when I is greater than I0. Ipar 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 byKd(PAR), is expressed as the form of (3) and calculated by using the SeaWiFS monthly mean data (averaged over 1998–2001) in January. Kd(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].
(3)Kd(PAR)=0.0864+0.884Kd(490)-0.00137[Kd(490)]-1,Kd(490)=0.0166+0.0773[chl]0.6715.


### 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 C (carbon) (mg m−3) through an empirical hyperbolic formula developed by Semovski and Wozniak [32]: C=ρmax(Chla/(Chla+K1/2))chla, where ρmax=90, and the half-saturation constant K1/2=0.477. Then a Redfield ratio (6.625) (Redfield et al. [33]) is adopted to convert C 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 forZ and D. 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 P, Z, and D decrease exponentially with depth in the initialization and take P as an example: P(k)=P(1)·e-(z(k)-z(1))/zch(k=1,2,…,14), where zch=90 m, which is characteristic scale depth.


## 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 global0.5°×0.5° 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(2)f(I,z)=II0·exp(1-II0),I=Ipar·exp(Kpar·z),
where I0 is the light intensity for maximum photosynthesis rate [30]. This function represents that with the increasing of I, photosynthesis rate increases exponentially when I is less than I0 and decreases rapidly when I is greater than I0. Ipar 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 byKd(PAR), is expressed as the form of (3) and calculated by using the SeaWiFS monthly mean data (averaged over 1998–2001) in January. Kd(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].
(3)Kd(PAR)=0.0864+0.884Kd(490)-0.00137[Kd(490)]-1,Kd(490)=0.0166+0.0773[chl]0.6715.


## 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 C (carbon) (mg m−3) through an empirical hyperbolic formula developed by Semovski and Wozniak [32]: C=ρmax(Chla/(Chla+K1/2))chla, where ρmax=90, and the half-saturation constant K1/2=0.477. Then a Redfield ratio (6.625) (Redfield et al. [33]) is adopted to convert C 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 forZ and D. 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 P, Z, and D decrease exponentially with depth in the initialization and take P as an example: P(k)=P(1)·e-(z(k)-z(1))/zch(k=1,2,…,14), where zch=90 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 byA-grid in horizontal direction and Z-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:(4)CF=12∑s,tWs(xs,t-x-s,t)2,
where x and x- are the modeled surface phytoplankton values and the observations, respectively. s refers to space while t refers to time. Ws is the spatial weight of the observation, which is defined as Ws=As/A-s, where As is the area of the grid, and A-s is the mean area of these grids which have observations.By introducing Lagrange multipliers, the adjoint model, equations of which are shown in AppendixB, 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 KPk+1=KPk+αkdk. In this paper, k is the assimilation step, dk is the direction (dk=(-V·5%/86400)·Gk→/|Gk→|, G→k is the gradient of cost function with respect to control variable), and αk 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)αk=1-0.1·(k-1),k=1,2,3,…,10,αk=0.1-0.01·(k-10),k=11,12,…,19,αk=0.01-0.001·(k-19),k=20,21,…,28.


## 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 R is selected to determine observations of how big an area (influencing zone) can be influenced by the parameter values of the independent grid. Eii,jj is the parameter of independent grid (ii, jj), and Fi,j is that of grid (i, j). Fi,j can be represented through linear interpolation of Eii,jj as the following relation [34]:
(6)Fi,j=∑ii,jjϕi,j,ii,jj·Eii,jj,ϕi,j,ii,jj=wi,j,ii,jj∑ii,jjwi,j,ii,jj,wi,j,ii,jj=R2-ri,j,ii,jj2R2+ri,j,ii,jj2,
where ϕi,j,ii,jj is the coefficient of linear interpolation, wi,j,ii,jj is the weighting factor, and ri,j,ii,jj is the distance between the independent grid (ii, jj) and the grid (i, j). After Eii,jj is tuned, Fi,j 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 valuesCV0 of the ecological parameters listed in Table 1, two types of parameter variations are constructed:
(7)Type1:CV(lat)=(-((lon-190)2+lat2)32400+1.25)·CV0,(8)Type2:CV(lat)=(((lon-190)2+lat2)32400+0.75)·CV0,
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.


### 3.1. Twin Experiment 1: Distribution Schemes of the Independent Grids

Distribution schemes of the independent grids are discussed by estimating spatially varyingVm in this section. Give the spatial variation to Vm 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 Vm.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, andR is 5°.(B)
There are 36 × 17 independent grids selected uniformly over the entire model domain, each from a 10° × 10° area, andR is 10°.(C)
There are 24 × 11 independent grids selected uniformly over the entire model domain, each from a 15° × 15° area, andR is 15°.Results of the three twin experiments above are shown in Table2. It is obvious that the more independent grids are used, the better the assimilation results are.Table 2
Comparison of the three strategies.
Scheme
Number of efficient independent grids
V
m RE (%)
MAE(mmol N m−3)
RCF
(A) 72 × 33
1731
0.21
0.0003
4 × 10−3
(B) 36 × 17
446
0.56
0.0006
6 × 10−3
(C) 24 × 11
184
1.16
0.0013
3 × 10−2


### 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.Vm, Dz, e, Dp, and Gm 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 Vm, Dz, e and spatial variation type 2 to Gm, Dp, 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 Figure2. 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 y=0.9x+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.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.


## 3.1. Twin Experiment 1: Distribution Schemes of the Independent Grids

Distribution schemes of the independent grids are discussed by estimating spatially varyingVm in this section. Give the spatial variation to Vm 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 Vm.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, andR is 5°.(B)
There are 36 × 17 independent grids selected uniformly over the entire model domain, each from a 10° × 10° area, andR is 10°.(C)
There are 24 × 11 independent grids selected uniformly over the entire model domain, each from a 15° × 15° area, andR is 15°.Results of the three twin experiments above are shown in Table2. It is obvious that the more independent grids are used, the better the assimilation results are.Table 2
Comparison of the three strategies.
Scheme
Number of efficient independent grids
V
m RE (%)
MAE(mmol N m−3)
RCF
(A) 72 × 33
1731
0.21
0.0003
4 × 10−3
(B) 36 × 17
446
0.56
0.0006
6 × 10−3
(C) 24 × 11
184
1.16
0.0013
3 × 10−2


## 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.Vm, Dz, e, Dp, and Gm 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 Vm, Dz, e and spatial variation type 2 to Gm, Dp, 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 Figure2. 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 y=0.9x+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.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.


## 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 Section2.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 3
Scope of assimilation area (noted by A).Figure 4
Cartoon for the setting of KP in the model. KPs in area A are optimized based on the conclusion in Section3.1, and the gray points represent the independent grids. KPs in area C are constant listed in Table 1. KPs in area B, the intermediate region, are obtained from area A and C by linear interpolation.Figure5 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.Assimilation results in Section4.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 CF1 and MAE1; after assimilation, the values of CF and MAE are denoted by CF28 and MAE28. The statistical area is A as shown in Figure 3.
(a)
(b)Results averaged over the 72 periods; (a) is the ratio of cost function to its initial value, and (b) is the MAE versus assimilation steps.
(a)
(b)


### 4.2. Comparative Experiments and Analysis

After assimilation, we obtain the temporal and spatial variation of each KP, denoted by KP(i,j,t), where i,j refer to space while t 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 KP(i,j,t), we average the spatial variation fields in time domain and then obtain one field of spatial variation KPS,KPSk(i,j)=∑tKPk(i,j,t), k=1,2,3,4,5. KPS1, KPS2, KPS3, KPS4, and KPS5 are the spatial variations of Vm, Dz, e, Dp, and Gm, which are shown in Figure 7.Spatial variation of KPs averaged in time domain. The values are scaled by the corresponding parameter values listed in Table1.
(a)
(b)
(c)
(d)
(e)(
2
) Temporally Varying KP. For each parameter of KP(i,j,t), 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, KPTk(t)=∑i,jKP(i,j,t), k=1,2,3,4,5. KPT1, KPT2, KPT3, KPT4, KPT5 are the spatial variation of Vm, Dz, e, Dp, Gm, 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, KPCk=∑i,jKPSk(i,j)=∑tKPTk(t), k=1,2,3,4,5. KPC1, KPC2, KPC3, KPC4, KPC5 are the constant Vm, Dz, e, Dp, and Gm, 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: KPSTk(i,j,t)=KPSk(i,j)·KPTk(t)/KPCk.Obviously, the estimated KP shows great spatial and temporal variations in Figures7 and 8, respectively. We find that Vm, Dz, and e decrease while Gm and Dp increase and vice versa. So we divide the five parameters into two groups: (1) Vm, Dz, and e and (2) Gm, Dp. Vm is the maximum uptake rate of nutrient by phytoplankton. When Vm increases, the growth of phytoplankton will be faster. Dz is the mortality rate of zooplankton. When Dz increases, the amount of zooplankton will reduce; on the contrary, the amount of phytoplankton will increase. e is the remineralization rate of detritus. When e 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, Gm is zooplankton maximum grazing rate, and Dp is the mortality rate of phytoplankton. When Gm decreases, the growth of phytoplankton becomes faster, and when Dp 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.Figure 8
Time series of KPs averaged in space. The values are scaled by the corresponding parameter values listed in Table1.Figure 9
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.The same conclusion can be drawn from results of correlation analysis in Tables3(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 Vm, Dz, and e in group (1) with the correlation coefficients reaching 0.93, and a positive correlation between Gm and Dp 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.Table 3
(a) Correlation coefficients between any two parameters of spatially varying KP. (b) Correlation coefficients between any two parameters of temporally varying KP.
(a)
e
D
z
D
p
G
m
V
m
0.58
0.72
−0.68
−0.65
e
0.91
−0.88
−0.92
D
z
−0.98
−0.99
D
p
0.99
(b)
e
D
z
D
p
G
m
V
m
0.76
0.92
−0.91
−0.91
e
0.93
−0.94
−0.94
D
z
−0.99
−0.99
D
p
0.99Run 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 Section4.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(i,j,t) 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.


## 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 Section2.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 3
Scope of assimilation area (noted by A).Figure 4
Cartoon for the setting of KP in the model. KPs in area A are optimized based on the conclusion in Section3.1, and the gray points represent the independent grids. KPs in area C are constant listed in Table 1. KPs in area B, the intermediate region, are obtained from area A and C by linear interpolation.Figure5 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.Assimilation results in Section4.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 CF1 and MAE1; after assimilation, the values of CF and MAE are denoted by CF28 and MAE28. The statistical area is A as shown in Figure 3.
(a)
(b)Results averaged over the 72 periods; (a) is the ratio of cost function to its initial value, and (b) is the MAE versus assimilation steps.
(a)
(b)


## 4.2. Comparative Experiments and Analysis

After assimilation, we obtain the temporal and spatial variation of each KP, denoted by KP(i,j,t), where i,j refer to space while t 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 KP(i,j,t), we average the spatial variation fields in time domain and then obtain one field of spatial variation KPS,KPSk(i,j)=∑tKPk(i,j,t), k=1,2,3,4,5. KPS1, KPS2, KPS3, KPS4, and KPS5 are the spatial variations of Vm, Dz, e, Dp, and Gm, which are shown in Figure 7.Spatial variation of KPs averaged in time domain. The values are scaled by the corresponding parameter values listed in Table1.
(a)
(b)
(c)
(d)
(e)(
2
) Temporally Varying KP. For each parameter of KP(i,j,t), 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, KPTk(t)=∑i,jKP(i,j,t), k=1,2,3,4,5. KPT1, KPT2, KPT3, KPT4, KPT5 are the spatial variation of Vm, Dz, e, Dp, Gm, 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, KPCk=∑i,jKPSk(i,j)=∑tKPTk(t), k=1,2,3,4,5. KPC1, KPC2, KPC3, KPC4, KPC5 are the constant Vm, Dz, e, Dp, and Gm, 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: KPSTk(i,j,t)=KPSk(i,j)·KPTk(t)/KPCk.Obviously, the estimated KP shows great spatial and temporal variations in Figures7 and 8, respectively. We find that Vm, Dz, and e decrease while Gm and Dp increase and vice versa. So we divide the five parameters into two groups: (1) Vm, Dz, and e and (2) Gm, Dp. Vm is the maximum uptake rate of nutrient by phytoplankton. When Vm increases, the growth of phytoplankton will be faster. Dz is the mortality rate of zooplankton. When Dz increases, the amount of zooplankton will reduce; on the contrary, the amount of phytoplankton will increase. e is the remineralization rate of detritus. When e 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, Gm is zooplankton maximum grazing rate, and Dp is the mortality rate of phytoplankton. When Gm decreases, the growth of phytoplankton becomes faster, and when Dp 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.Figure 8
Time series of KPs averaged in space. The values are scaled by the corresponding parameter values listed in Table1.Figure 9
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.The same conclusion can be drawn from results of correlation analysis in Tables3(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 Vm, Dz, and e in group (1) with the correlation coefficients reaching 0.93, and a positive correlation between Gm and Dp 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.Table 3
(a) Correlation coefficients between any two parameters of spatially varying KP. (b) Correlation coefficients between any two parameters of temporally varying KP.
(a)
e
D
z
D
p
G
m
V
m
0.58
0.72
−0.68
−0.65
e
0.91
−0.88
−0.92
D
z
−0.98
−0.99
D
p
0.99
(b)
e
D
z
D
p
G
m
V
m
0.76
0.92
−0.91
−0.91
e
0.93
−0.94
−0.94
D
z
−0.99
−0.99
D
p
0.99Run 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 Section4.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(i,j,t) 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.



---

*Source: 373540-2013-05-27.xml*