Interactive comment on “ Reducing the model-data misfit in a marine ecosystem model using periodic parameters and Linear Quadratic Optimal Control ”

[11pt,a4paper]article amsfonts amsmath graphicx wrapfig graphicx amsmath graphicx epstopdf pdfpages tabularx algorithmic longtable lscape amsthm caption epstopdf multirow multicol url tabularx font=small,labelfont=bf,singlelinecheck=false float subfig color ruled plain equationsection Comments to referee 2: We thank you for your comments, and are sorry for the inconsistencies. We give answers here, referring to your numbering.


Introduction
In this paper we present an application of the LQOC (Linear Quadratic Optimal Control) method on a parameter optimization in a marine ecosystem model.This method (see for example Kwakernaak and Sivan, 1972;Lewis and Syrmos, 1995) is a mathematical technique to compute time dependent parameters (in some applications also called controls) in linear dynamical systems.The key goal of the method is to minimize a given quadratic cost function subject to a lin-ear system.For the cost function, a model-data fit in leastsquares formulation can be used, and thus the method is applicable to parameter optimization or model calibration problems.The optimal time-dependent parameters are obtained via an algorithm that basically uses the system matrices of the underlying linear model.When applying the method to a non-linear dynamical system (as most marine ecosystem models), the linearization of the system is essential.For this purpose, a reference trajectory of model variables and parameters is needed, which can be based on observational data and parameter guesses.By an appropriate choice of the parameter trajectory, periodic parameters can be obtained.These can then be used to improve the original non-linear model.
Marine ecosystem models describe biogeochemical processes in the ocean and are used, e.g., for calculating the effect of marine photosynthesis on the global carbon cycle.Typically, such kind of models have several parameters, for example, growth and mortality rates for the different types of plankton taken into account.Since most of these parameters are not known exactly and are difficult to measure, parameter identification or estimation is an important tool to calibrate a model.Parameter identification or estimation is usually done by performing a parameter optimization in order to minimize the misfit between model output and given data, commonly represented by a least-squares type cost functional.Additionally, uncertainty estimates corresponding to data errors may be computed.A parameter optimization may improve a model's quality also to another extent: If a model still shows deficiencies after the parameters have been optimized, reasons other than inappropriate parameter values are likely to be responsible for remaining poor model behavior, see for example Fasham and Evans (1995), Hurtt and Armstromg (1996), Fennel et al. (2001), Prunet et al. (1996).
The computational effort to perform such kind of optimization runs for the coupled system of ocean circulation and marine biogeochemistry that describes a marine ecosystem is quite high, especially in three space dimensions.Thus, several simplifications may be used: One of them is to compute the marine biogeochemistry in a so-called offline mode, i.e., to solve the transport equations for the tracers with precomputed ocean circulation fields (velocity, temperature, salinity) as forcing or input.Another approach is to use onedimensional models, which simulate a single water column only.This simplification is motivated by the fact that most of the ecosystem processes (as for example growth and dying) are happening locally in space and that the main spatial interactions are vertical mixing and sinking of organic matter.Moreover, it has been shown that parameters obtained in a one-dimensional optimization can also be beneficial when used in three-dimensional computations, see Oschlies and Schartau (2005).
The motivation for this paper is based on the results obtained for a typical marine ecosystem model, namely an NPZD model (N for dissolved inorganic nitrogen, P for phytoplankton, Z for zooplankton and D for detritus) introduced in Oschlies and Garc ¸on (1999).As was reported in several publications with different optimization algorithms, the quality of the model fit to observations was not optimal, and in some cases it was difficult to identify the parameters uniquely, see for example Ward (2009), Ward et al. (2010), Rückelt et al. (2010b).In most cases (and in all these studies), the parameters of the marine ecosystem models are assumed to be temporally constant.This reflects the aim to obtain a model that is applicable for arbitrary time intervals.In contrast, in our work we allow the parameters to vary temporally over the year while remaining periodic over all years of the considered time interval.
Our main research question is if such kind of relaxation is able to significantly improve the model-to-data fit.Eknes and Evensen (2002) and Schartau et al. (2001) have examined the possibility of using a sequential data assimilation method for state estimation in a biological model.On the other hand, there are several papers on parameter estimation only, see Schartau et al. (2001), Fasham and Evans (1995), Hurtt and Armstromg (1996), Fennel et al. (2001), Prunet et al. (1996), Matear (1995), Spitz et al. (1998).Work by Losa et al. (2003) combined state and parameter estimation using a sequential weak constraint parameter estimation in an ecosystem model.An example for time-dependent parameters is introduced in the work by Mattern (2012), where a statistical emulator technique to estimate time-dependent values for two parameters of a 3-dimensional biological ocean model is used.The author demonstrated that emulator techniques are valuable tools for data assimilation and for analyzing and improving biological ocean models.He also used temporally changing parameters, but without imposing annual periodicity.
We use the annual periodicity constraint on the parameters in order to allow for some temporal flexibility of the parameters but at the same time to retain the temporal universality of the optimized model, e.g., for application to time periods outside the range of observations.The LQOC used here does not require a prescribed periodic parameterization of the parameters, it will automatically generate an optimal periodic function for each parameter.Moreover, it allows to balance the two aims of a good model-data fit on one hand and parameter periodicity on the other by introducing weighting matrices.
We verify our approach in validation and prediction experiments employing the optimized periodic parameters in the original non-linear model.
The structure of the paper is as follows: in Sect. 2 we briefly describe the model, the data and the cost function which we use for optimization.In Sect. 3 we briefly describe the LQOC method and its application on the NPZD model.In Sect. 4 we present our results obtained by the LQOC method with respect to the quality of the model-data misfit and the periodicity of the parameters.We furthermore show results for validation and prediction with the original nonlinear model using the optimized periodic parameters.Section 5 ends the paper with some conclusions.

Model equation and optimization problem
The model used here, as example, is a one-dimensional marine ecosystem model presented in Oschlies and Garc ¸on (1999).It is of NPZD type, i.e., it simulates the interaction of dissolved inorganic nitrogen (N), phytoplankton (P), zooplankton (Z) and detritus (D), whose concentrations (in mmol N m −3 ) are denoted by the model variables (y l ) l=N, P, Z, D =: y.These four variables are functions y l : [0, t e ]×[−H, 0] → R of space and time, with H denoting the depth of the water column and t e the total integration time.
The model is then given as the following system of partial differential equations: here z denotes the vertical spatial coordinate.The q l are the biogeochemical coupling terms which depend on space and time via light intensity and also on temperature, and on most of the parameters summarized in the vector u.
In this spatially one-dimensional setting, the only physical process taken into account is vertical diffusion, which appears as a space and time-dependent mixing coefficient κ, taken (as well as temperature) from the Ocean Circulation and Climate Advanced Model OCCAM (see Sinha and Yool, 2006) in hourly profiles.The equation for detritus also contains a sinking term with constant speed ω D > 0, which is also optimized as u 12 , whereas ω N = ω P = ω Z = 0 in Eq. ( 1).In total, we have p = 12 parameters u = (u 1 , . . ., u p ) to be optimized.The biogeochemical coupling (or source-minus-sink) terms are taken from Oschlies and Garc ¸on (1999) and read (2) Here J is the daily averaged phytoplankton growth rate as a function of depth z and time t, and G is the grazing function: u 7 +u 6 y P 2 .The circulation data (taken from an ocean model) are the turbulent mixing coefficient κ = κ(z, t) and the temperature T = T (z, t), which is used in the non-linear term c T where c = 1.066 is kept constant.Table 1 lists the model parameters with their original symbols as in Oschlies and Garc ¸on (1999).For more details see also Schartau and Oschlies (2003a).

Observational data and corresponding model output
The observational data used here, denoted by y obs , is taken from the Bermuda Atlantic Time-series Study (called BATS) as a part of the US Joint Global Ocean Flux Study, see Michaels and Knap (1996).
The BATS data are provided by the Bermuda Biological Station for Research (BBSR) situated in the Atlantic Ocean, 700 miles from the East Coast of the US at coordinates 31 • N 64 • W (see http://bats.bios.edu/).
In this work, each observational data has to be compared to an equivalent value generated by the model.For this purpose, the model output is interpolated in time and space to match the observational data.In addition, some transformations have to be done: modeled zooplankton has to be integrated in space, and chlorophyll a values are calculated by multiplying the nitrogen-based concentration of the modeled phytoplankton by a factor of 1.59 mg Chl/(mmol N), which corresponds to a chlorophyll to carbon mass ratio of 1 : 50 and a C : N mole ratio of 106 : 16.For more details see Ward (2009).
Summarizing, there are five types of observational data y obs := (y l,obs ) l=N, P, Z, D, PP , which correspond to aggregated values ȳ := ( ȳl ) l=N, P, Z, D, PP of the model output.The used data and their corresponding model variables are 1.Dissolved inorganic nitrogen y N, obs (mmol m −3 ) corresponding to model variable y N =: ȳN .
Here using a constant conversion factor of 1.59 mg (Chl a)/(mmol N).
5. Carbon fixation or primary production (PP) as carbon uptake y PP, obs in mmol C m −3 d −1 .As modeled primary production, the temporal mean of the model output y P multiplied by the phytoplankton growth rate J (y N )y P .

The optimization problem
The aim of the optimization is to fit the model output ȳ that was aggregated, in the above mentioned way, to the given observational data y obs over a chosen time interval of j max years.We denote by N l,j the number of observational data for y l,obs for each observed quantity l = N, P, Z, D, PP in year j = 1, . . ., j max .These numbers may be different for each l and j .The i-th observational in year j of y l,obs is denoted by y l,obs j,i , and the corresponding aggregated model output value by ( ȳl j,i ).We now firstly compute the averaged annual misfit per model output/tracer, weighted using the inverse of the standard deviations taken from Schartau and Oschlies (2003a) and summarized in the vector σ = (σ l ) l=N, P, Z, D, PP = (0.1, 0.01, 0.01, 0.0357, 0.025), and by the number N l,j of observational per tracer and year: , l = N, P, Z, D, PP, j = 1, . . ., j max .
If there are no observational data for a state variable/tracer in a year (i.e., N l,j = 0), the sum is set to zero.The overall cost function is then calculated as where N total is the total number of non-zero terms F l,j actually occurring in the sum.In the usual case we have N total = 5j max .If ever N l,j = 0 and thus F l,j = 0 for a year and tracer, N total is decreased accordingly.
Table 1.Parameters of the ecosystem model to be optimized with the LQOC method.Here u 0 = (u 0,i ) i=1,...,12 is the vector of parameters taken from Oschlies and Garc ¸on (1999), min(u i ) and max(u i ) their respective upper and lower bounds used in Schartau and Oschlies (2003a). Parameter

Application of linear quadratic optimal control to the NPZD model
In this section we apply the LQOC method to the discretized version of the NPZD model.The LQOC method is widely used in engineering applications and well studied from the mathematical side, see e.g., Anderson and Moore (1971), Casti (1987), Lunze (1997), Sima (1996).Non-linear problems can be treated by linearization, see Clemens (1993).We present the details of the linearization and enforcement of the periodicity of the parameters.

Discretization scheme
We here give a brief description of the temporal and spatial discretization of the model equations described in Eq. ( 1).
The vertical grid consists of 32 layers with thickness increasing with depth.In Rückelt et al. (2010a) it has already been demonstrated that, at least from the point of view of the optimization results, the vertical model grid can be reduced to this number instead of the originally employed 66.It has been demonstrated that optimization of the model yield practically identical results w.r.t.parameter match and quality of the optimal solution.Anyhow, the method used here can as well be applied for the original 66 layers.The model is discretized in time using an operator splitting method: given a time-step size t (one hour in the model), the discretized scheme reads (5) Here is the vector of all four tracers and u k the parameter vector (which is here already assumed to be time-varying), both at the current time step k.The total number of time steps is M. The matrices L diff k , L sink k are 4×4 block-diagonal and represent the discretization of diffusion (discretized by second order central differences) and sinking (discretized by an upstream scheme), respectively, and I is the identity matrix.
The interpretation of the scheme (Eq.5) is the following: In every time step k → k + 1, at first the non-linear coupling terms q k = (q N k , q P k , q Z k , q D k ) are computed at every spatial grid point and integrated by four explicit Euler steps with step size t 4 , each of which is described by the operator Then, an explicit Euler step with full step-size t is formed for the sinking term, represented by the matrix [I + tL sink k ].This matrix does only depend on the time step k if the sinking velocity w s is to be optimized.Finally, an implicit Euler step is applied for the diffusion operator, discretized with second order central differences.Due to κ = κ(z, t) the resulting matrix [I − tL diff k ] depends on the current time step k.The discrete system can now be formally written as where f is a non-linear function.

Reference tracer trajectory
In order to apply the LQOC method to the discretized nonlinear system (Eq.7), we perform a linearization around reference trajectories of state y ref and control u ref .The latter is described in Sect.3.5.The relation between model variables and observational data described in Sect.2.1 can be used to define the reference trajectory y ref as follows 1. y N, ref := y N, obs .
3. y Z, ref is taken as constant over the whole water column such that the integral equals the value of the observational data y Z, obs (which is a mean value).Here the formula (Eq. 3) is used.
Because the data for the different tracers are not given for the same instances of time, we only use those instances where data for all four tracers are available.Only for those instances reference values for all tracers can be computed by the above transformations.To describe this procedure, we define for every year j = 1, . . ., j max : -N l,j : number of data for y l,obs , l = N, P, Z, D.
t l j,i : instance of time of y l,obs j,i , i = 1, . . ., N l,j , rounded to an integer number.
-I j : set of instances where data for all tracers are available: -N j : number of elements in this intersection.
Because we have no reference value y ref at the end of each year, we use linear interpolation between the data y obs j,N j and y obs j +1,1 to generate such value, denoted by y ref j,N j +1 .Since there are no data for every vertical layer, we additionally interpolate the data linearly in space to obtain a reference trajectory in every grid point.

Linearization
We apply the linearization on each time interval [t j,i , t j,i+1 ] separately.For this purpose, we split the interval into sections of the length t (here t = 1 h), as shown in Fig. 1.
The discrete time steps τ k in this interval are They correspond to the steps used in the temporal discretization, see Sect. 3. Now we linearize the model around the data y ref j,i+1 and the parameter vector u ref k (described below in Sect.3.5).For this purpose, we introduce QOC method to the discretized nonerform a linearization around refer-  discrete time steps τ k in this interval are They correspond to the steps used in the temporal discretization, see Section 3.Here x k ∈ R m is the deviation of the model output (i.e., all tracers over the whole water column) from the reference trajectory, v k ∈ R p is the deviation of the actual parameters from the reference parameter trajectory.
The dimension m of x k is determined by the number of tracers (in our case 4) and the numbers of grid cells in the water column (in our case 32), which results in m = 4 • 32.The dimension p of v k is determined by the number of model parameters, here p = 12.The matrices A k ∈ R m×m and B k ∈ R m×p are called system matrix and input matrix, respectively.For this example, they were generated by Algorithmic or Automatic Differentiation (AD), see Griewank (2000).Here we used the software TAF (Transformations of Algorithms in Fortran), see Giering and Kaminski (1998).
Summarizing, the state equation reads This is the typical form of a discrete linear system with state variable x k and parameter (or control) v k .

Application of the LQOC theory
The theory of linear quadratic optimal control gives a formula for the optimal parameter trajectory v that minimizes the cost function under the constraint Eq. ( 8).Here for every k -Q k ∈ R m×m is a positive semi-definite diagonal weighting matrix for the state vector, -R k ∈ R p×p is a positive definite diagonal weighting matrix for the parameter vector.
The matrices Q k and R k are usually chosen to be diagonal.They reflect the relative importance of keeping tracer variables and parameters, respectively, close to their reference trajectories.In our case, this translates to quality of the model-to-data fit and periodicity of the parameters, see the next subsection.Selecting large elements of either matrix will emphasize the corresponding effect.
At end of the optimization on all intervals [t j,i , t j,i+1 ], j = 1, . . ., j max , i = 1, . . ., N j , we obtain the optimal state variable and the optimal periodic parameters (y * k , u * k ), for k = 1, . . ., M − 1.The realization of the LQOC method is presented in algorithm 1.

Choice of the reference parameter trajectory
A main objective of this work is to enforce periodicity of the parameters.To obtain annually periodic parameters, we denote the length of a time period -which in our case is one year -measured in time steps by T = t 1,N 1 +1 .We now choose the reference trajectory for the parameter to be Here u 0 ∈ R p is an initial guess, in our case we took the values from Oschlies and Garc ¸on (1999), compare Table 1.These values are used in the first year (k ≤ T ) only.As a result, in the first year the choice of the cost function (Eq.9) will force periodicity to the constant reference parameters.This effect can be reduced by choosing appropriate small values in the matrices R k in the first year.In the following years, the difference of the current parameter u k to its counterpart from the year before is minimized.This enforces periodicity.
Thus the crucial point in adjusting the matrices R k throughout an optimization run is to allow both: -for sufficiently large deviation from the constant reference parameters in the first year, to enable their temporal variation, -and for smaller deviation and thus more or less strict periodicity in the following years.

Particular choice of Q k , R k
In our example, the weighting matrices Q k are taken as constant for all k, namely where σ l are the standard deviations taken from the original cost function (Eq.4).The matrices R k are taken as These values are chosen differently in the first year (on one hand) and in all subsequent years (on the other hand).In all years except the first one (i.e., k ≥ T ), the R k are used to enforce periodicity of the parameters.The bigger the r k n for these years are, the better periodicity of the parameters is expected.Following this idea, the optimal choice for the R k in the first year would be just zero matrices.But, by this choice, the requirements for the LQOC method and algorithm 1where the R k have to be positive definite -are not satisfied.As a consequence, it is desirable to chose the r k n for the first year as small as possible.
On the other hand, the choice of the r k n in the first year can be used to keep the parameters in the admissible bounds: Since they can be forced to keep in the vicinity of the initial guess in the first year and to stay close due to the periodicity enforcement in the following years, by a careful setting of the r k n in the first year both aims (periodicity and boundedness) can be balanced.This effect is not guaranteed by the LQOC method, but turned out to be realizable in our case, see the next section.
Summarizing, for our computations we chose for n = 1, . . ., p, the values where u 0 is as listed in Table 1.

Optimization results
In this section we present the results of the parameter optimization runs performed with the LQOC method.We show both the obtained fit of the linearized model output (with optimal parameters) to the data and the annual periodicity of the obtained parameters.Note that we only compare values of the original cost function F , the function J actually minimized in the LQOC setting is just a tool of the method.

Fit of linear model output to observational data
This section shows a comparison between the optimized model output, obtained by the LQOC method with periodic parameters, and the observational data.As a reference we also compare the results to those obtained by a direct optimization of the original non-linear model using constant parameters with a Sequential Quadratic Programming (SQP) method that takes into account parameter bounds.This method was used in Rückelt et al. (2010b).
Figure 2 shows the model results, obtained by the LQOC method with periodic parameters, for aggregated model output ȳ and the observational data y obs for the years 1994 to 1998 for the uppermost layer at depths z ≈ 5 m.Shown are results for a part of the whole time interval at some distinct depth layers only.The total number of depth layers considered in the optimization process is 32 and the total number of time steps is 43 800.In contrast to the results obtained for constant model parameters with the original non-linear model, the LQOC method with periodic parameters gives a nearly perfect fit of the data.Even substantial concentration changes that occur between some neighboring observational data points (e.g., for y N, mod , in 1994, 1995 or 1997) can be captured by the optimized trajectory.We performed the Algorithm 1. Algorithm of the LQOC method.for j = 1, . . ., j max do for i = 1, . . ., N j do (a) Given the equation system as , and the cost function as 1. Solve the first matrix difference Riccati equation with final condition P t j,i+1 = Q t j,i+1 .

Solve the second matrix difference Riccati equation
with final condition h t j,i+1 = 0.
3. Solve the two auxiliary matrix difference equations 5. Obtain the optimal parameters vector u end end optimization for the years 1994 to 1998, in contrast to the years 1991 to 1996 that were used in Rückelt et al. (2010b), since no zooplankton data are available at BATS for the years 1991 to 1993.This would be disadvantageous for the linearization procedure in the LQOC method.
In Rückelt et al. (2010b), a minimum value of the cost function (Eq.4) of F ≈ 70 was obtained for optimized constant parameters for the time 1991 to 1996.For the time 1994 to 1998, the value obtained by the SQP method and constant parameters is very similar.Also the quality of the fit -depicted in Fig. 2 -is comparable.The better fit obtained by the LQOC method also results in a significantly lower value (F ≈ 1.35) of the original cost function (Eq.4). Figure 3 shows the mismatch between model output and the reference data -which are interpolated values of the sparse observational data -for all points in time and space.
Figure 3 also shows -except for dissolved inorganic nitrogen ȳN -a better fit at the surface than in deeper layer.The possible reason for this is the lack of observational data on the lower layers, which requires interpolation over relatively large space intervals.Thus, the interpolation error there is bigger than in upper layers where the database is denser.This, naturally, affects the quality of the reference trajectory, and thus also of the linearized model and the obtained optimal parameters.
Figure 4 shows the LQOC results with periodic parameters and the observational data y obs for the years 1994 to 1998 for the lower layer at depths z ≈ 184.32 m.Here the LQOC method with the periodic parameters provides a nearly perfect fit, in contrast to the one obtained with constant parameters.

Sensitivity with respect to the weighting matrices R k
To examine the effect of the weighting matrices R k in the first year, on the behavior of the parameters and the cost function F , we have additionally performed sensitivity experiments with different entries r k i of the weighting matrices R k for k ≤ T .We present two additional experiments for n = 1, . . ., p with The values of min(u n ) and max(u n ) are listed in Table 1.The corresponding values of the r k n for these two choices and the r k n from Eq. ( 11) (here called reference simulation) are shown    in Table 2.The trajectories of the tracers y, the parameters u, and the value of the cost function F , depend heavily on the choice of the corresponding entries r k n in the matrix R k .Figure 5 shows the trajectories for three tracers and different r k n .All experiments show only minor differences from the reference simulation with the r k n from Eq. ( 11).The results show that a decrease in the entry r k n can lead to a small decrease in the cost function.The sensitivities of the parameters with respect to the choice of r k n can be seen in Fig. 6.It is obvious that for smaller values of r k n the variability of the parameters is getting larger.

Periodicity of the parameters
In this section we show that the above model-to-data fit can be achieved with parameters that are almost annually periodic.Enforcement of periodicity was achieved by an appropriate adjustment of the matrices R k in the cost function (Eq.9) used in the LQOC framework, see Sect. 3. It was also possible to keep the parameters in their desired bounds (see Table 1), although the LQOC method does not to impose these bounds explicitly.
Figures 6 illustrates the temporal behavior of the parameters that were optimized with the LQOC method.Depicted are only the ten that show a temporal variation.Two parameters remain constant in time.These figures show different trajectories for each parameter for two years with the different choices of the r k n , compare Table 2.As mentioned above,  it is obvious that for a smaller r k n , the amplitude of the parameters increases, but it always remains almost periodic.Since the periodicity of the parameters is nearly perfect, then it is enough to plot for 2 yr.
The parameters controlling growth of phytoplankton, namely the maximum growth rate a and the initial slope of the P-I curve α, show in Fig. 6 both maximum values in early summer and in winter, with a clear minimum value  n , the reference simulation (dashed), with larger r k n (solid) and with smaller r k n (dotted) for the upper layer at depth z ≈ 5 m and years 1994-1998.in spring during the peak of the annual chlorophyll signal.This is consistent with earlier assimilation studies that, for assumed constant parameters, tended to overestimate plankton production at BATS during the bloom end of winter and, at the same time, tended to underestimate production in oligotrophic summer conditions and in early winter, see Schartau et al. (2001).Such a trend to relatively high values of α has also been found in earlier studies that optimized parameter values by data assimilation, see Fasham and Evans (1995), Schartau et al. (2001).Earlier studies assuming timeindependent parameter values have attributed relatively high values of α to the absence of a dial cycle in the turbulent mixing, which might allow for substantial phytoplankton growth even in winter during reduced daytime mixing, see Schartau and Oschlies (2003a).This is consistent with the findings of the current study, that also suggest high values of α during the period of deep mixing in winter.In addition, our optimized model predicts even higher values of the initial slope parameters α for late spring and early summer, where the mixed layer is usually shallow and growth is limited by nutrients rather than light in the surface mixed layer.A large value of α can, however, help to establish a subsurface chlorophyll maximum in better agreement with the observations.This was also noted by Schartau and Oschlies (2003b).Our results reported here indicate that high values of α may, at BATS, be more important for the establishment of the deep chlorophyll in late spring than for the maintenance of phytoplankton production during periods of deep mixing in winter.Maintenance of high primary production during summer has been difficult to achieve by earlier models run at BATS (Schartau et al., 2001).As nutrient supply to the surface waters is low during the stratified season, models with fixed carbon-to-nutrient stoichiometry and constant model parameters do not seem to be able to reach observed levels of primary production in the surface layer, see Schartau and Oschlies (2003b).In the current study, the carbonto-nutrient factor used to convert simulated (nitrogen-based) primary production to observed (carbon-based) primary production is constant as well.However, the seasonally varying parameters can contribute to maintain high levels of primary production during summer in the absence of substantial inputs of new nutrients.This is realized by enhanced recycling of biomass, evident by high maximum grazing rates, high assimilation efficiencies, high prey capture efficiencies and high zooplankton excretion in late spring and early summer.Similarly, remineralization of detritus is highest in late spring as well.These high rates all contribute to fast recycling of nutrients in the surface ocean, which helps to maintain observed high rates of primary production and thereby reduces the model-data misfit function.The relative deviations of the calculated parameters are shown in Fig. 7.

Validation of the non-linear model with periodic parameters
In this section, the periodic parameters obtained by the LQOC method are used in a validation experiment using the original non-linear NPZD model.We run the original nonlinear model using these parameters without further optimization for the years 1994 to 1998 and analyze the corresponding model-data misfit.
Figure 8 shows a comparison of the model output using optimal periodic parameters and optimal constant parameters (obtained by the SQP method), as well as the observational data in the uppermost layer.The use of periodic parameters in comparison to constant parameters results in a significantly better model-data fit.The results for y P and y D are almost perfect, whereas for the other two tracers they are slightly worse.The results look similar for all layers.
In order to allow for a quantitative comparison between our results and those obtained with constant parameters by Rückelt et al. (2010b), we give the corresponding values of the original cost function (Eq.4) in Table 3.The better fit with periodic parameters also results in a significantly reduced value F ≈ 15.05 of the cost function (Eq.4).We are not able to obtain the same value of the cost function as with the linearized model, see Sect.4.1, which is reasonable since there we have used the observational data for the reference trajectory.
The quality of the obtained periodic parameters depends on the length of the time intervals [t j,i , t j,i+1 ] on which the linearization (with a constant reference value for the model output) is applied.If no data are available as reference y ref k , these intervals have to be enlarged, which presumably will reduce the quality of the optimized parameters.This in turn will influence the quality of the validation.Fig. 6.Periodicity of optimal parameters (uk,n)n=1,...,10 obtained by the LQOC method, min(u),max(u) are respective the upper and lower bounds listed in table 1, u opt is the optimal parameter obtained by Rückelt et al. (2010).We point out that min(u) and max(u) are not shown in all plots.Because the values are very low or high.Fig. 6.Periodicity of optimal parameters (u k,n ) n=1,...,10 obtained by the LQOC method, min(u), max(u) are, respectively, the upper and lower bounds listed in Table 1, u opt is the optimal parameter obtained by Rückelt et al. (2010b).We point out that min(u) and max(u) are not shown in all plots.Because the values are very low or high.Fig. 8. Observational BATS data y l,obs ,l = N,P,Z,D and aggregated model trajectories ȳl ,l = N,P,Z,D, obtained by using the optimal periodic parameters in the non-linear model, and the optimal trajectories for the constant parameters with a Sequential Quadratic Programming (SQP) method at depth z ≈ 5m for the years 1994 -1998.
Fig. 8. Observational BATS data y l,obs , l N, P, Z, D, and aggregated model trajectories ȳl , l = N, P, Z, D, obtained by using the optimal periodic parameters in the non-linear model, and the optimal trajectories for the constant parameters with a sequential quadratic programming (SQP) method at depth z ≈ 5 m for the years 1994-1998.

Prediction experiment
In this section we present a prediction experiment with the optimal periodic parameters in the original non-linear NPZD model.We now use only a part (two years, 1994 and 1995) of the time interval to determine optimal periodic parameters using the LQOC method.Then we use these parameters on the remaining part (three years, 1996 to 1998) of the time interval.In these three years, the periodic parameters obtained during the first two years are applied without further optimization.
Figure 9 shows a comparison between the predicted model output and the observational data for the years 1996 to 1998 in the uppermost layer.The fit is quite good.The qualitative behavior of the tracers at different times and spatial layers is similar.Table 3 shows the resulting values of the cost function.
The fit of the predicted output is slightly worse than for the output in the validation experiment described in Sect.4.4, but still much better than the results obtained with the optimized constant parameters.

Conclusions
In this paper, we use the method of linear quadratic optimal control (LQOC) to determine optimal periodic parameters in a one-dimensional marine ecosystem model.We demonstrate that the LQOC method can be applied on the considered parameter optimization problem for a non-linear NPZD type model using a linearization technique around reference  trajectories of model variables (biogeochemical tracers) and parameters, where the system matrices were obtained by automatic differentiation.
We show how a reference tracer trajectory can be built from sparse data, and how an appropriate choice of reference parameter trajectory can be used to obtain annual periodic parameters that additionally stay in prescribed bounds.
The linearized model obtained in this way gives a very good model-data fit with almost perfect annually periodic parameters.Even with the available small number of observational data typical to oceanographic time series sites, the quality of the fit is very high.Specifically, it is much better than the one previously obtained by optimization of the nonlinear model with fixed model parameters.
The obtained periodic parameters are used in the original non-linear model.Using them, the model is able to reproduce and predict the real data much better than the non-linear model with optimized constant parameters.
The method allows to further analyze temporal deviations of individual parameters about the annual mean.This may help in making inferences about processes that the model cannot describe well when constant parameters are used.This latter analysis should contribute to a better understanding of model deficiencies and, eventually, help to improve marine ecosystem models.

Fig. 1 .
Fig. 1.Example for an interval on which the linearization is applied.

Fig. 1 .
Fig. 1.Example for an interval on which the linearization is applied.

Fig. 2 .
Fig. 2. Observational data y l,obs ,l = N,P,Z,D and aggregated model trajectories ȳl ,l = N,P,Z,D , optimized with periodic parameters obtained by the LQOC method and with a Sequential Quadratic Programming (SQP)method.Values are shown for the upper layer at depth z ≈ 5m and years 1994-1998.

Fig. 2 .Fig. 3 .
Fig. 2. Observational data y l,obs , l = N, P, Z, D, and aggregated model trajectories ȳl , l = P, Z, D, optimized with periodic parameters obtained by the LQOC method and with a sequential quadratic programming (SQP)method.Values are shown for the upper layer at depth z ≈ 5 m and years 1994-1998.

Fig. 3 .Fig. 4 .
Fig. 3. Model-to-data misfit for four tracers with respect to space and time obtained by the LQOC method.

Fig. 4 .
Fig. 4. Observational data y l,obs , l = N, P, Z, D, and aggregated model trajectories ȳl , l = N, P, Z, D, optimized with periodic parameters obtained by the LQOC method.Values are shown for the lower layer at depths z ≈ 184.32 m.

Fig. 5 .
Fig.5.Model output trajectories with different r k n, the reference simulation (dashed), with larger r k n (solid) and with smaller r k n (dotted) for the upper layer at depth z ≈ 5m and years1994-1998.

Fig. 5 .
Fig. 5. Model output trajectories with different r kn , the reference simulation (dashed), with larger r k n (solid) and with smaller r k n (dotted) for the upper layer at depth z ≈ 5 m and years1994-1998.

Fig. 7 .
Fig. 7. Relative temporal variations of some of the model parameters in a typical year.Since the periodicity of the parameters is nearly perfect, no difference between the five years is visible.

Fig. 7 .
Fig. 7.Relative temporal variations of some of the model parameters in a typical year.Since the periodicity of the parameters is nearly perfect, no difference between the five years is visible.

Fig. 9 .
Fig.9.Observational BATS data y l,obs ,l = N,P,Z,D and aggregated model trajectories ȳl ,l = N,P,Z,D, obtained by using the optimal periodic parameters in the non-linear model as a prediction at depth z ≈ 5m for the years 1996 -1998.

Fig. 9 .
Fig. 9.Observational BATS data y l,obs , l = N, P, Z, D, and aggregated model trajectories ȳl , l = N, P, Z, D, obtained by using the optimal periodic parameters in the non-linear model as a prediction at depth z ≈ 5 m for the years 1996-1998.
ref and control u ref .The latter is y Z,ref − y P,obs .separately.For this purpose, we split the interval into sections of the length ∆t (here ∆t = 1h), as shown in Figure1.The

Table 2 .
Values of r n i for the reference simulation (second column) and the two additional sensitivity experiments.

Table 3 .
Values of the cost function for each year: for the optimization (with constant and periodic parameters), the validation presented in Sect.4.4, and the prediction presented in Sect.4.5.