Journal of Marine Systems 45 (2004) 1 – 20 www.elsevier.com/locate/jmarsys

Weak constraint parameter estimation for a simple ocean ecosystem model: what can we learn about the model and data? Svetlana N. Losa a,b,c, Gennady A. Kivman a,b,*, Vladimir A. Ryabchenko b b

a Alfred-Wegener-Institut fu¨r Polar and Meeresforschung, Postfach 120161, 27515 Bremerhaven, Germany St. Petersburg Branch of P.P. Shirshov Institute of Oceanology, Russian Academy of Science, 30, Pervaya Liniya, 199053 St. Petersburg, Russia c Department of Oceanography, Dalhousie University, Halifax, Nova Scotia, Canada B3H 4J1

Received 27 September 2001; accepted 25 August 2003

Abstract Weak constraint parameter estimation is applied to tune a 0-D 4-component (phytoplankton, zooplankton, detritus, dissolved inorganic nitrogen) model for the ecosystem of the upper mixed layer in the North Atlantic. The model is constrained by the monthly mean Nimbus-7 chlorophyll data. The basin under consideration is separated into 5
5j boxes where the model is optimized independently. Adjusted model parameters exhibit significant spatial variations. We also show how to evaluate a posteriori the relative credibility of the data and the model. A method for estimating the accuracy of model parameterizations is suggested. For Bermuda Station ‘‘S’’, the major contribution to model errors comes from uncertainties in parameterizations for the phytoplankton mortality, zooplankton losses and breakdown of detritus. D 2003 Elsevier B.V. All rights reserved. Keywords: Parameter estimation; Ecosystem modeling; Data assimilation

1. Introduction A common feature of ocean ecosystem models in existence is that all involve a big number of poorly known parameters which are difficult, if impossible, to measure directly. Thus, the problem of tuning these parameters and calibrating models of the ocean ecosystems becomes of great importance. Variational data assimilation (VAR) provides a systematic tool to * Corresponding author. Alfred-Wegener-Institut fu¨r Polar and Meeresforschung, Postfach 120161, 27515 Bremerhaven, Germany. Fax: +49-471-4831-1797. E-mail address: [email protected] (G.A. Kivman). 0924-7963/$ - see front matter D 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.jmarsys.2003.08.005

tackle this problem. At the heart of the method is the search for values of the model parameters which minimize a distance (least-square, as a rule) between the system trajectory and data. In doing so, the model equations may be imposed either as strong or as weak constraints. To the best of the authors’ knowledge, all of the numerous attempts to optimize ecosystem model parameters (Fasham and Evans, 1995; Matear, 1995; Prunet et al., 1996; Hurtt and Amstromg, 1996; Spitz et al., 1998; Fennel et al., 2001; Schartau et al., 2001) share two common features. First, data were assimilated into local (0-D and 1-D) models and thus the horizontal advective transport of biological substances

2

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

was ignored. Second, the model equations were treated as the strong constraints. These raise two important questions: – are these estimates applicable at other locations? – are the model errors truly negligible? There is some evidence for spatial and temporal variability of the ecosystem model parameters (Platt et al., 1991). In a recent paper, Hemmings et al. (2003) showed ‘‘that in some cases calibrations obtained for one domain can be successfully applied to another, but that this is not true generally.’’ On the other side, the neglected model uncertainties are interpreted in the course of strong constraint parameter estimation as model parameter variations and may result either in unrealistic parameter estimates (Matear, 1995; Hurtt and Amstromg, 1996) or in a solution deviating far from the data (Fasham and Evans, 1995; Fennel et al., 2001). Remaining in the context of 0-D models, our goals are to investigate the spatial variability of the model parameters and model errors. To this end, we employed weak constraint parameter estimation. It is worth noting that, as a rule, weak constraint data assimilation has been used for state estimation. However, if the model parameters are poorly known, state estimation with fixed model parameters can produce unacceptable results (Kagan et al., 1997). Weak constraint parameter estimation is a combination of the adjoint method and the generalized inversion, thus, it allows one to take into account errors inherent in the model and data as well as find optimal values of the poorly known model parameters. Previously this method has been applied to linear or weakly nonlinear systems (Brummelhuis et al., 1993; Eknes and Evensen, 1997; Gong et al., 1998). Natvik et al. (2001) mentioned that a combined parameter and state estimation problem could be formulated for strongly nonlinear ocean ecosystem models. This is the data assimilation technique we utilize. Certainly, for the study aimed at an investigation of spatial variability of model parameters, it might be better to perform the optimization in a 2D or 3D context when accounting for the advection of biological components, but this generalization dramatically increases the dimension of the control space. Instead, we consider an ensemble of local box models covering the North Atlantic with 5j resolution and

believe that weak constraint data assimilation takes into account the errors caused by neglecting the transport of the biological components. Though this assumption still remains rather restrictive, it is much weaker than that underlying strong constraint data assimilation. As was emphasized by Fennel et al. (2001), optimization of the model parameters can play an important role in model development and verification and allows one to test the model’s skill to reproduce known data sets. The superiority of the weak constraint method to the strong constraint approach in recovering a system trajectory that much better fits the data has been shown in several studies (see, for example, (Bennett and Thorburn, 1992; Gong et al., 1998; Kivman et al., 2001). We will demonstrate that in addition, weak constraint data assimilation makes it possible to derive valuable supplementary information about the model itself. This information is completely lost when the strong constraint scheme is used. Namely, we present a method of making a conclusion of whether the observational information or the model is closer to the truth. It is worth while to emphasize that the inverse solution obtained with weak constraint data assimilation does not satisfy the model equations exactly, and thus, it is difficult to infer the correct biogeochemical fluxes between different model compartments, which is a particular goal of ecosystem modelling and data assimilation (Evans, 1999). In this paper, we propose a secondary inversion procedure that allows us to redistribute model equation residuals between different biogeochemical processes and, given the system trajectory that is close to the data but does not satisfy the model equations exactly, to obtain estimates of the fluxes. What is probably even more important is that this procedure offers a tool to detect the major source of the model errors and thus determine the path to improving model performance. The paper is organized in the following manner. The data assimilation procedure is considered in the next section. Section 3 describes the ecosystem model for the upper mixed layer (UML). Experiment design and results for the North Atlantic are described in Sections 4 and 5. Results of an identical twin experiment performed to illustrate controllability of the system by chlorophyl data are presented in Section 6. Section 7 is devoted to an evaluation of the

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

credibility of the model equations, the accuracy of the data, and the quality of the parameterizations of the involved biological processes. Section 7 contains a summary and discussion.

2. Weak constrained parameter estimation Let us write the ecosystem model in the form dc ¼ Mp ðcÞ; dt

ð1Þ

where the vector c contains unknown concentrations of model components, Mp is a nonlinear model operator, and p is the set of model parameters. Following Bennett (1992), we present observational information about the ecosystem as Lm ðcÞ ¼ dm ;

m ¼ 1;: : :; M ;

ð2Þ

where Lm (c) are observational functionals acting on the system state variables c, dm are the data, and M is amount of observational information. The generalized inversion (Bennett, 1992) consists in minimizing the cost function: Fðc; pÞ ¼ x

Z T 0

þb

2 dc  Mp ðcÞ dt dt

M X

ALm ðcÞ  dm A2 ;

ð3Þ

m¼1

where (dc/dt  Mp(c)) are residuals referred to the model equations, (Lm(c)  dm) are data misfits, x and b are some weights that measure the reliability of the model and data, respectively. Since our goal is to study the climatological dynamics of the ecosystem, the periodicity condition cð0Þ ¼ cðT Þ;

ð4Þ

is imposed. Here T is a 1-year period. The cost function (Eq. (3)) is equal to the log-likelihood of the probability density function (PDF) qðc; pAd1 ;: : :; dM Þ ¼ A expðFðc; pÞÞ

ð5Þ

3

Here A is a normalization constant. Under the assumption that the data and model errors are Gaussian with zero mean and variances of (2b) 1/2 and (2x) 1/2 (Tarantola, 1987), Eq. (5) expresses the state of our knowledge about the model variables and parameters when the data d1,. . .,dM are observed. Thus, the minimizer of the cost function gives the most probable realization of the system trajectory and the model parameters. Generally, the errors in the equations for different model compartments may correlate (they also may correlate in time). In this case, instead of a single weight x in Eq. (3), a non-diagonal matrix will appear. Here we adopt the simplest hypothesis: the errors in different model equations and different time increments are independent and identically distributed. The solution minimizing the cost function approximately satisfies the model equations and fits the data. It is worth noting that treating the model inputs p as variables of the cost function (Eq. (3)) in addition to the system state variables c is a distinguishing feature of the approach used. The strong constraint parameter estimation commonly used in ecosystem data assimilation is a limit case of the formulation presented when x ! l, and thus, the cost function depends only on p. Imposing dynamics as a weak constraint has a very important implication. Namely, the cost function becomes an analytical function of all its variables p and c, that is not the case for the strong constraint data assimilation. As a result, the cost function varies more regularly over the control space and, as was demonstrated by Evensen and Fabrio (1997), the minimization problem becomes mathematically better posed compared to the strong constraint formulation. Choosing proper values for the weights x and b is a very important issue which strongly affects the quality of the solution. Small b yield a system trajectory diverging from the data. On the other hand, as our experience has shown, assigning large values to b, that is equivalent to overtrusting the data, can lead to negative values of optimized model parameters and/or some model variables. In principle, the weights should be chosen to be inversely proportional to the squared variances of the model and the data errors, however, the model uncertainties can be hardly esti-

4

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

mated a priori. In general, we know little about the model uncertainties, or even about the observational errors. This is especially true for the satellite-derived surface chlorophyll concentrations assimilated in this study. It is only known that the relative errors may be in the range of 15– 35% for the open ocean (Balch et al., 1992), however, they may exhibit notable spatial variations and rise in highly clouded areas (Antoine et al., 1996) and in the coastal regions where high concentrations of the yellow substance are often observed (Evans and Gordon, 1994). Therefore, it is rather difficult, if impossible, to choose the data error variance a priori. Instead, we adopt estimating the weighting a posteriori. Currently, there are a number of approaches to estimating the weights a posteriori: the maximum likelihood method (Cramer, 1954), the generalized cross-validation (Wahba, 1990), etc. Here we utilize the maximum data cost criterion (MDC) recently put forward by Kivman et al. (2001). From a computational point of view, very large or small weights will lead to ill-conditioning of the problem. To explain the idea underlying the MDC, let us consider the second term (referred to as the data cost) of the cost function (Eq. (3)). If we compute its value for an inverse solution c*(b) that corresponds to some b and treat it as a function of b, it can be shown that the data cost tends to be 0 as b ! 0 and as b ! l, and here the data cost should have a maximum. The MDC suggests choosing the value of b which corresponds to the maximum of the data cost. Our experiments with different b revealed that choosing values significantly higher than that recommended by the MDC resulted in an inadmissible inverse solution (either negative model variables or model parameters). As the minimizer of the cost function depends only on the b/x ratio, hereinafter x is taken as 1. Thus, the problem to be solved takes the form

ðc; pÞ ¼ Arg min

(Z  T 0

þb

M X

2 dc  Mp ðcÞ dt dt ) 2

ALm ðcÞ  dm A

ð6Þ

After discretizing equations (Eq. (1)) in time, we obtain the cost function

Fðc; pÞ ¼

with b chosen according to the MDC.

Dt

k¼0

þb

M X

i2  Mp ðck Þ Dt

ALm ðcÞ  dm A2 Þ;

ð7Þ

m¼1

defined over a finite dimensional space, where c = {c0,. . ., cK}, and K is the number of time steps per 1 year. Following Evensen and Fabrio (1997) and Kagan et al. (1997), we apply the conjugate gradient descent method to minimize the cost function (Eq. (7)) directly. Given a current iteration for the control variables and the gradient of the cost function which is computed analytically, the next iteration is readily determined without forward and backward model integrations.

3. Ecosystem model The model is based on the four-component model of Popova et al. (1997) and includes phytoplankton P, zooplankton Z, detritus D and dissolved inorganic nitrogen N as dependent variables. We write the model equations describing the processes of entrainment – detrainment at the lower UML boundary, sinking of detritus with velocity wg, and inner biological sources and sinks Bi as h

dci þ qih0 þ di3 wg c3  hBi ¼ 0; dt

ð8Þ

where h is the UML thickness, c1 = P, c2 = Z, c3 = D, c4 = N, qi (i = 1,2,3,4) are the turbulent fluxes at the lower UML boundary, and di3 is the Kroneker symbol. The fluxes qi are parameterized in terms of the entrainment velocity (we=dh/dt) in a common way, i.e.

;

m¼1

K h X ckþ1  ck

qih0 ¼

8 < we ðci  c*i Þ; if we > 0; :

0;

if we V0;

ð9Þ

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

where ci* are the component concentrations at the upper boundary of the seasonal pycnocline. Except for c4* ¼ a1 þ a2 h;

ð10Þ

for some coefficients a1 and a2 of a linear vertical profile of the dissolved inorganic nitrogen in the seasonal pycnocline, all other ci* were specified equal to 0. The biological sources and sinks Bi are expressed as 8 if i ¼ 1; PP  G1  De1 > > > > > > > > if i ¼ 2; < b1 G1 þ b2 G2  De2 Bi ¼ > > > ð1  b1 ÞG1  b2 G2  De3 þ De1 if i ¼ 3; > > > > > : PP þ mDe2 þ De3 if i ¼ 4; ð11Þ where PP = JQP is the daily mean phytoplankton growth rate, J is the light-limited growth rate, Q is a non-dimensional limiting factor, G1 and G2 are grazing rates of the zooplankton on the phytoplankton and on the detritus, respectively, b1 and b2 are equivalent assimilation efficiencies, De1 is the rate of phytoplankton natural mortality, De2 is the rate of zooplankton losses due to excretion, natural mortality and grazing of higher order predators, De3 is the rate of detritus breakdown, and m is a fraction of the total zooplankton losses transformed into dissolved inorganic nitrogen. The remaining part 1  m of zooplankton losses transforms into detritus which is considered to be instantly exported to the seasonal pycnocline. Following Fasham (1993), J, Q, De1, De2 and De3 are given as J¼ Q¼

2 h

Z sZ 0

N k1 þN

h

F½I0 ðtÞexpfðkx þ kc PÞzgdzdt;

0

; De1 ¼

l1 P 2 k5 þP

; De2 ¼

l2 Z 2 k6 þZ

; De3 ¼ l4 D; ð12Þ

where F(I) = VpaI/(V p2 + a2I 2)1/2, Vp = V *p 1.066Tw is the maximum phytoplankton growth rate, a is the initial slope of the P – I curve, I = I(z,t) is the photosynthetically available radiation at depth z, I0 = I(0,t) is assumed to be proportional to the absorbed total

5

irradiance at the sea surface, and the ratio kp of I0 to the irradiance being prescribed; t is measured in days and is 0 at sunrise and s at noon; kw is the light attenuation coefficient due to water, kc is the phytoplankton self-shading parameter, Tw is the water temperature, z is the vertical coordinate, k1, k5, k6, l1, l2, l3 and Vp* are model parameters. The zooplankton grazing rates G1 and G2 are described by the following expressions: gr1 ZP2 ; k3 þ r1 P2 þ r2 D2 gr2 ZD2 G2 ¼ ; k3 þ r1 P2 þ r2 D2

G1 ¼

ð13Þ

which are multi-prey generalizations of the Holling type III function. Here r1 and r2 are the weighted preferences for phytoplankton and detritus, g and k3 are the maximum specific grazing rate and the halfsaturation constant for grazing.

4. Experiment design We have applied the assimilation method for the ecosystem model to simulate the seasonal cycle of the UML ecosystem and to estimate the poorly known model parameters in each cell of the 5
5j grid covering the North Atlantic from 30jN to 60jN. In doing so, we minimize the cost function (Eq. (6)) with simultaneous estimation of the data weight b independently for each grid cell. The monthly mean solar irradiance at the ocean surface was taken from Oberhuber (1988). The UML depth was estimated from temperature profiles (Levitus and Boyer, 1994) as the depth at which the temperature is 0.5jC less than that at the surface. The coefficients a1 and a2 in Eq. (10) were calculated as the linear regression coefficients for the vertical nitrate profile from nitrate data at the standard depths (Conkright et al., 1994). The data constraining the model were the monthly mean values of the surface chlorophyll concentration observed from the satellite Nimbus-7 by the coastal zone color scanner (CZCS) and averaged over the period from 1979 to 1985. All the external forcing data and chlorophyll data were averaged over each 5
5j box.

6

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

To convert the chlorophyll concentration that is not involved in the model to the phytoplankton biomass, the chlorophyll to nitrogen ratio, Chl/N (or the chlorophyll to carbon ratio Chl/C) for phytoplankton was used. The ratio is spatially and temporally variable and depends on numerous factors (Geider et al., 1996; Hurtt andArmstrong, 1996). To account for some of them, the following expression (Cloern et al., 1995) Chl=C ¼ a þ be

cT x dI

e

N ; k1 þ N

Table 1 Initial model parameters Symbol

Parameter

l1

phytoplankton maximum specific mortality rate half-saturation constant for nutrient uptake phytoplankton mortality half-saturation constant initial slope of the P – I curve phytoplankton maximum growth rate constant light attenuation coefficient by water absorption phytoplankton self-shading coefficient ratio of photosynthetically available radiation to total irradiance zooplankton maximum loss rate maximum specific grazing rate half-saturation constant for grazing zooplankton loss rate half-saturation constant zooplankton assimilation efficiencies zooplankton feeding preferences detrital breakdown rate detrital sinking velocity nitrogen fraction of zooplankton losses

k1 k5 a V *p

ð14Þ

was used. Here a, b, c, and d are some constants. Prescribing an averaged value of the C/N ratio ( = 6.625), we immediately come to the expression for Chl/N. This dependence of the chlorophyll concentration on the model variables ( P,N) and forcing (I,Tw) was used to calculate Lm(c) in Eq. (2).To account for the fact that the signal measured by the satellite reflects information of the chlorophyll within a thin upper layer, I is attenuated over the upper 35 m (Vinogradov et al., 1995) in Eq. (14). A crucial issue is the number of the model parameters which could be estimated from the 12 monthly mean chlorophyll data. Since we do not penalize the deviations of the parameters from the first guess values, clearly we cannot tune more than 12 parameters. Even this number may be too high and result in a very poorly constrained problem. Instead, we decided to tune only six parameters and keep the others fixed. Popova (1995) studied the sensitivity of the model trajectories to parameter variations and found the sensitivity (measured as the relative root-mean square deviation between the reference solution and that computed after adding a zero mean Gaussian noise with the relative standard deviation of 50% to reference values of the model parameters) to be quite different for different model compartments and different locations. In particular, the phytoplankton, zooplankton and detritus are much more sensitive to the model parameters variations than dissolved inorganic nitrogen. For the Bermuda Station ‘‘S’’ and the Ocean Weather Station ‘‘I’’ which are located in our region of interest, these model compartments are most sensitive to g, l1, l2, a and V *, p k3, respectively. Consequently, these six parameters (the maximum phytoplankton growth rate constant V *, p the phyto-

kw kc kp

l2 g k3 k6 b1, b2 r1, r2 l4 wg m

Values

Units

0.05

day 1

0.3

Mmol N m 3

0.2

Mmol N m 3

0.025

m2 W 1 day 1

0.6

day 1

0.04

m 1

0.03

m2 Mmol N 1

0.41

0.25

day 1

0.73

day 1

0.1

Mmol N m 3

0.2

Mmol N m 3

0.75 0.7,0.3 0.05 10 0.8

day 1 m day 1

plankton maximum specific mortality rate l1, the zooplankton maximum loss rate l2, the zooplankton maximum ingestion g, the zooplankton ingestion half-saturation constant k3, and the initial slope of the P – I curve a) can be derived from the data with the highest accuracy, and so there were the parameters we adjusted. The values of Popova (1995) summarized in Table 1 were used as the first guess.

5. Results Figs. 1 and 2 depict respectively the horizontal distributions of the annual mean and temporal evolution of zonal mean chlorophyll a concentration in the UML obtained with and without data assimilation

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

7

Fig. 1. Annual horizontal distribution of the chlorophyll a concentration (mg m3) in the UML of the North Atlantic: (a) the model solution, (b) the inverse solution and (c) the CZCS satellite data averaged over 1979 – 1985.

against the satellite data. The maps were obtained by taking all the quantities at the centers of the boxes and plotting isolines, using linear interpolation. As it is seen, the model solution fits the data much better in the southern part of the basin than it does in the mid-latitudes that is in agreement with Popova (1995).She found the model to be more adequate for the Bermuda Station ‘‘S’’ (32jN, 64jW) compared to the Ocean Weather Station ‘‘I’’ (59jN, 19jW) and suggested that data assimilation could markedly improve model’s skill in the mid-latitudes. Indeed, the inverse solution is capable of reproducing a northward increase of the chlorophyll concentration, with its local maximum around 52jN, 32jW, and a significant, but not insufficient, increase near the European

coast (Fig. 1). Spring and autumn blooms (timing and magnitude) in the mid-latitudes as well as the phenomena of ‘‘rolling carpet’’ (a lag in timing of the spring bloom with increasing latitude) are also distinctly present in the inverse solution (Fig. 2). However, the inverse solution still misses an increase in chlorophyll concentration in winter in the mid-latitudes and underestimates the autumn bloom in the high latitudes. The latter discrepancies may, at least, partly be explained by the acquisition of the ocean surface color data: it was not systematic and some regions were persistently cloudy. Using interpolation procedures to fill missing information can lead to overestimating chlorophyll concentrations at the middle and high latitudes of the ocean, starting from 40jN

8

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

Fig. 2. The time evolution of zonal mean chlorophyll a concentration (mg m3) in the UML of the North Atlantic. Panels (a), (b), (c) correspond the same solutions and data as in Fig. 1.

(Antoine et al., 1996). We will show below that the data assimilation method used is able to discern the quality of these data. The annual mean distribution of dissolved inorganic nitrogen (Fig. 3) is in close agreement with the nitrate data from the World Ocean Database (2001) and reproduces the northward increase as well as a wavy structure of the dissolved inorganic nitrogen spatial distribution in mid-latitudes between 40jN and 50jN. It is worth noting that the

effect of data assimilation on dissolved inorganic nitrogen is marginal, which is in agreement with Popova (1995) who found a very low sensitivity of this model compartment to parameter variations. Fig. 4 depicts the spatial distributions of the mean annual maximum phytoplankton growth rate Vp = Tw V *1.066 and optimized model parameters a, g, p k3, l1, l2. To compare our results with estimates of the phytoplankton growth rate obtained from obser-

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

9

Fig. 3. Annual horizontal distribution of dissolved inorganic nitrogen (Mmol N m3) obtained from the inverse solution for the UML of the North Atlantic.

vational data by Platt et al. (1991, referred to as PCS91), we present in Fig. 2a the phytoplankton growth rate Vp rather than the parameter Vp* which was optimized. As is seen from the comparison of Fig. 4a,b with the mean annual values of Vp and a obtained by averaging seasonal estimates of PCS91 Table 2, our results are in qualitative agreement with these experimental estimates. Namely, in both cases, Vp and a increase as we move southward. However, our values of Vp exceed the estimates of PCS91 everywhere. The optimal values of a are overestimated in the southern part and underestimated in the northern part of the North Atlantic, while being in perfect agreement with the PCS91 estimate for the central part of the basin. A point that should be emphasized is that we present here a comparison of our estimates, averaged over the subdomains, with those presented in PCS91, averaged over the annual cycle. One can notice rather high spatial variances of our estimates within each subdomain. The spatial mean of the parameters is not necessarily the best estimate from the point of view of model performance. The same is true for the annual mean. According to PCS91, Vp and a can vary 2- and 10-fold (respectively) during the annual cycle. Hence, we can also expect large standard deviations of the estimates presented in Table 2. In any case, the estimates of PCS91 are within the confidence intervals of our estimates.

We can also compare the optimized values obtained for a and Vp with the more recent estimates of Kyewalyanga et al. (1998, referred to as K98) (see Table 3) for subdomains I– V depicted in Fig. 4., and we see that Vp is in very close agreement with K98. The estimates of a are somewhat overestimated, but remain rather close to the spring values of K98, except in the region III. Spatial distributions of the other adjusted model parameters (Fig. 4c – f) are also noticeably varied, however, the variability stays within ranges adopted in the modelling community (Table 4). One can see some prevalence of the parameter variations in the meridianal direction, in particular, the zooplankton parameters g, k3 and l2 increase and the phytoplankton maximum specific mortality rate l1 decreases in the northward direction. In addition, the zooplankton parameters increase in the northeastern corner of the region. The spatial distribution of the parameters is, obviously, a result of a combined effect of several factors such as solar irradiance, temperature, etc., which may affect latitudinal changes of chemical conditions and plankton species composition. However, it is rather difficult to distinguish which of the physical, chemical and biological factors, and in which region, contributes more to the spatial variability of the physiological model parameters. Kyewalyanga et al. (1998) found some correlations between the a and Vp and so-called environmental and biological variables,

10

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

Fig. 4. Horizontal distributions of optimized model parameters: (a) the maximum phytoplankton growth rate Vp, mg C (mg Chl h) 1; (b) the initial slope of the P – I curve a, mg C m2 (mg Chl h W) 1; (c) the maximum specific rate g, day 1, and (d) half-saturation constant k3, mmol m 3, for grazing; the maximum specific mortality rate (day 1) for phytoplankton, l1, (e) and zooplankton, l2, (f). For Vp and a regions where independent estimates of K98 are available are plotted.

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

11

Fig. 4 (continued).

however, these correlations varied significantly and unsystematically from one biogeochemical province to another. The significant spatial variations of the adjusted model parameters point out that it is unlikely to find a unified set of values of the model parameters for the North Atlantic, or even more for the global ocean.

6. Twin experiment It was mentioned above that the minimizer of the cost function (Eq. (3)) gives the most probable realizations of the system trajectory together with a set of

Table 2 Estimates of the maximum phytoplankton growth rate Vp, and the initial slope of the P – I curve a Latitudes

2

mg C m a,mg Chl Wh

Vp,

mg C mg Chlh

PCS91 Current estimates PCS91 Current estimates 51 – 70jN 0.082 38 – 50jN 0.063 11 – 37jN 0.094

0.043 F 0.023 0.066 F 0.040 0.258 F 0.185

2.4 3.0 4.0

4.7 F 3.8 4.5 F 2.0 7.2 F 3.4

the model parameters conditioned by the data available. The results presented in the previous section express the most we can say about the ecosystem and revealed regular spatial variability of the optimized model parameters having only the CZCS data to work with. If the data assimilated could indicate too little about the system, the minimum of the cost function would be very flat. In this case, the minimization would resemble a random walk covering a rather large area around the optimal values of the model parameters. Since the minimization is performed in each grid box independently, the spatial structure of the estimated model parameters would be very irregular that is not the case. However, we still have to admit a possibility that though the inverse solution is closer to the data better than the initial model trajectory, this better fit to observations might be achieved at the expense of corrupting some unobserved variables (either model compartments or some parameters) of the cost function. Therefore, we have to investigate the quality of the inverse solution. The problem of how successful data assimilation can be if only the chlorophyll data are assimilated has

12

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

Table 3 Values of the maximum phytoplankton growth rate Vp, and the initial slope of the P – I curve a for a few biogeochemical provinces from K98 Biogechemical provinces

Parameter values 2

mg C m a,mg Chl Wh

I II III IV V

Vp,

mg C mg Chlh

Spring

Autumn

Spring

Autumn

0.073 F 0.048 0.100 F 0.045 0.075 F 0.028 0.078 F 0.025 0.069 F 0.032

0.019 F 0.007 0.040 F 0.020 0.037 F 0.006 0.049 F 0.037 0.023 F 0.008

3.30 F 2.63 6.01 F 2.35 6.88 F 3.30 8.23 F 2.62 6.64 F 3.37

2.24 F 1.56 4.51 F 2.09 4.58 F 2.32 4.43 F 2.68 4.88 F 1.28

Hessian depends on many factors such as the data values, the first guess parameters, etc., and may behave uniquely in different regions of the control space. Thus, no general conclusion can be drawn and, as is pointed out by Schartau et al. (2001), an identical twin experiment is better to be performed after, not before obtaining the inverse solution constrained by the real data. This ensures that the twin experiment is not conducted within a completely different region of the control space and thus gives us an insight on the error bars of the inverse solution. For the twin experiment, we picked up the box 50– 55jN, 35 –40jW and took the inverse solution and the corresponding set of the model parameters as the truth. To generate synthetic data, we added Gaussian noise with a zero mean and a relative standard deviation of 30% to the true values of the monthly mean chlorophyll concentrations computed for the reference solution. Then, we repeated data assimilation, following the scheme described above and starting from the same first guess. Table 5 summarizes the results of the parameter estimation. It is seen that we are able to improve the estimates for all parameters, although the errors for k3 and g still remain relatively high. A source of these errors may be that both parameters are involved in the computation of the same term (the grazing rate) and, even knowing it

two facets in the context of our study: how well can we recover the whole system state and how good are the parameter estimates? The first problem was addressed by Natvik et al. (2001) and Eknes and Evensen (2002). Natvik et al. (2001) reported that they were able to recover a reference solution with a weak constraint inverse ecosystem model when the distance between measurements was 5 days. When this distance was set to 17 days, the conjugate gradient converged to a local minimum of the cost function. An analogous result was reported by Eknes and Evensen (2002) for the Ensemble Kalman filter. Srokosz (1998) studied the accuracy of the parameter estimation and compared different strong constraint methods for assimilation of the chlorophyll data. He also concluded that it was possible to recover the model parameters with high accuracy from the chlorophyll data if the data were not too noisy. Formally, the full answer to these questions is given by analysis of the Hessian of the cost function (Fennel et al., 2001) which provides the formal error bars for the inverse solution. However, the gradientbased minimization does not output the Hessian directly, and therefore, the accuracy of the inverse solution is investigated in most studies by identical twin experiments. The problem is that since ecosystem parameter estimation problems are nonlinear, the Table 4 Ranges of parameter values Parameter

G k3 l1 l2

Our estimates

Other estimates

min

max

mean

min

0.27 0.02 0.01 0.09

1.42 1.30 0.09 0.32

0.71 F 0.24 0.28 F 0.36 0.04 F 0.03 0.22 F 0.04

0.30 0.06 0.015 0.075

max (JGOFS Rep., 1997) (Popova et al., 1997) (Prunet et al., 1996a,b) (JGOFS Rep., 1997)

2.06 1.38 0.10 0.36

(Spitz et al., 2001) (Prunet et al., 1996a,b) (JGOFS Rep., 1997) (Evans and Parslow, 1985)

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20 Table 5 Twin experiment Parameters

l2

g

Vp

l1

k3

a

First guess True value Optimised value

0.25 0.23 0.23

0.73 0.83 0.75

0.6 1.2 1.0

0.05 0.007 0.007

0.1 0.9 1.3

0.025 0.007 0.006

precisely, one cannot determine these two parameters independently. The model trajectory markedly underestimates the chlorophyll (Fig. 5) and it is no surprise that assimilation of the chlorophyll data yields an inverse solution much closer to the truth than the pure model solution. The lack of the chlorophyll is partly corrected by decreasing the phytoplankton specific mortality rate l1 by a factor of about 7 (see Table 5). The

13

improvement of the system trajectory is more pronounced for the second half of the year where the synthetic data appeared to be of better quality and almost coincide with the truth. Assimilating the chlorophyll data also allowed us to correct the deficit of the dissolved inorganic nitrogen exhibited by the pure model, solution, another factor resulting in higher chlorophyll concentrations. The major deviations of the inverse solution from the truth, though quite small, occur between April and June when the assimilated data are less accurate. The detritus computed without data assimilation stays rather close to the reference solution except for the period between February and April. Constraining the solution with the chlorophyll data results in a much better fit to the reference solution for February and March. The error in April

Fig. 5. Identical twin experiment for the 50 – 55j N 35 – 40j W box. The dashed curve is the initial solution; the solid curve is the reference solution; the circles are monthly mean reference chlorophyll; the stars are synthetic chlorophyll data; the thick solid curve is the inverse solution.

14

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

still remains rather large and is caused by inaccuracy of the April mean chlorophyll concentration. The model trajectory for the zooplankton is similar to the reference solution (especially for the winter and late autumn) and assimilation of the relatively inaccurate data for April – June does not improve the fit to the true system trajectory for this period. However, the inverse solution becomes almost indistinguishable with the true zooplankton concentrations for July – September when the very precise data come into play. We also repeated the data assimilation, keeping fixed the model parameters optimised for the box located between 30 – 35jN and 60 – 65jW, i.e. the cost function (Eq. (7)) depended only on the state variables in this experiment. With this setup, the inverse solution had little in common with the reference system trajectory(we do not present it here). This result implies that errors in specification of the model parameters cannot be compensated by accounting for the model noise in data assimilation which is in agreement with Kagan et al. (1997). General conclusions which can be drawn from the twin experiment performed are the following. First, assimilating only the chlorophyll data improved the estimate for all model parameters and compartments compared to the first guess values. Second, increasing the data quality has significant positive impact on the fitness of the inverse solution to the truth and this effect is localized in time, i.e. the better the observations for a particular period, the closer the inverse solution is to the truth. Third, deviations of the optimized values of the model parameters from the truth obtained in the identical twin experiment turned out to be much smaller than the spatial variability of the model parameters presented in the previous section. Thus, this spatial variability is a stable feature and not a random response to changing the forcing fields and data due to a lack of constraints, to get a sharply piked PDF (Eq. (5)).

7. Inference about the model and the data Traditionally, data assimilation is viewed as a tool for fitting models to data, however, when the model equations are imposed as weak constraints, we obtain something more. Computing the model equation residuals, we can glean additional valuable informa-

tion about the model. In particular, we will show below how it is possible to infer something about the relative quality of the model and the data and about the major source of the model uncertainties. 7.1. Which is better: the model or the data? In the simplest case when we have two observations y1 and y2 of the same quantity x, the best linear unbiased estimate x˜ of x is given by the minimizer of the cost function IðxÞ ¼ ðx  y1 Þ2 þ bðx  y2 Þ2 ;

ð15Þ

where b is equal to the ratio of the squared error variances of the observations y1 and y2 (Ghil and Malanotte-Rizzoli, 1991). Clearly ð˜x  y1 Þ2 bð˜x  y2 Þ2

¼ b;

ð16Þ

holds for the minimizer x˜ of I(x). Such a simple equality cannot be derived in the multidimensional nonlinear case considered here, but still we can use the ratio of the terms in the right-hand side of Eq. (6) P b M ALm ðcÞ  dm A2 r ¼ R h m¼1 i2 T dc 0 dt  Mp ðcÞ dt

ð17Þ

to characterize the relative errors in the data and model solution. This non-dimensional quantity is more suitable for the analysis than b in Eq. (6) which has the dimension of m 2d. The spatial distribution of r (Fig. 6) depicts that the satellite data are less accurate than the model (r < 1) in the coastal and shelf regions, especially in the northeastern part of the North Atlantic. This feature can be easily explained. The yellow substance and colored sediments or detritus presented in these regions strongly distort the signal coming from the ocean surface to the satellite, due to absorption and dispersion, and can lead to systematic overestimation of the chlorophyll concentration (Evans and Gordon, 1994). Moreover, the north-eastern part of the North Atlantic is known to be highly clouded. This last circumstance is one of the reasons for overestimation of the satel-

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

15

Fig. 6. Horizontal distribution of the ratio r of the data and model residuals.

lite-derived chlorophyll concentrations in this region (Antoine et al., 1996). On the other hand, the model equations cannot be considered to be perfect. As Fig. 6 depicts, r < 10 for most of the North Atlantic and we cannot neglect the contribution of the model noise to the cost function. It would be very useful to determine the main sources of the model errors which may contain a systematic component (bias). One can expect worse model performance in regions of strong ocean currents since the advection of the biological substances is neglected. Fig. 7 presents the annual mean model equation residuals normalized by the corresponding biological sources: residuals of the first model equation are normalized by primary production, PP; residuals of the second equation are normalized by zooplankton production, b1G1 + b2G2; residuals of the third equation are normalized by De1 phytoplankton mortality; residuals of the fourth equation are normalized by the m-fraction of zooplankton losses and detritus breakdown to dissolved organic nitrogen, mDe2 + De3. As can be seen, all the equations involve systematic errors which are markedly higher in the regions of the Gulf Stream and the North Atlantic current. Thus, the four-compartment ecosystem model we employ properly describes the dynamics of the ecosystem only in the central part of the ocean, far from the strong boundary currents. In other regions, additional factors and processes controlling behavior of

the biogeochemical system have to be taken into consideration. Generally speaking, advective transport of biological substances contributes mainly to systematic errors in reality, especially in regions of jet currents. This may not be the case in our study where data are averaged over 5
5j boxes, however. Under such averaging effects, advection may more weakly manifest and other model restrictions such as the timeinvariant linear dissolved inorganic nitrogen profile and absence of zooplankton in the seasonal pycnocline may play an important role. Moreover, to improve the performance of the model it has to be augmented with additional ecosystem components affecting the seasonal dynamics of the phytoplankton and the ecosystem as a whole. For example, other groups of phytoplankton and zooplankton, dissolved organic matter, and bacteria are not included in the model. Finally, some the model parameterizations used can contribute significantly to the model errors. 7.2. Inference about the model parameterizations and fluxes A particular goal of ecosystem modelling and data assimilation is to determine biogeochemical fluxes (Evans, 1999). Under ideal conditions, we would ˜ , N˜) which covers all possess a data set c˜ = (P˜, Z˜, D model state variables at each time increment. How could we estimate fluxes in this case? The most straightforward option is to use a particular ecosystem

16

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

Fig. 7. Annual model equation residuals normalized by the total biological source in the equations for phytoplankton (a), zooplankton (b), detritus (c) and dissolved inorganic nitrogen (d).

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

model, based on some parameterizations of the fluxes, to optimize the model parameters involved and to find those values resulting in the model trajectory. This will give us the best model trajectory in the sense of fitting the data. We can then compute the fluxes from this solution by means of explicit mathematical formulas used in the model to parameterize biological processes, however, if this trajectory does not fit the data within the data error (which is generally the case), we cannot claim that the model derived fluxes are close to the truth when the model trajectory is not. Another feasible scheme is to compute the fluxes directly from the data by means of the same mathematical formulas, but with the optimized model parameters. One will immediately find, however, that the temporal evolution of the biological compartments, as given by the data, is in contradiction with the fluxes estimated by these means. At this point, we have to realize that we are dealing with an inverse problem of estimating the fluxes from two different sorts of information (the model and the indirect data), neither of which is perfect. Thus, we can prescribe them as weak constraints, however, there is a general property of mass balance which is more fundamental than any parameterization and we wish to preserve it. In the case of the model considered, if we denote the data covering the whole ˜ , N˜), the following assimilation period as c˜ = (P˜, Z˜, D equations must be satisfied PP  G1  De1 ¼

dP˜ þ h1 q1 ; dt

b1 G1 þ b2 G2  De2 ¼

dZ˜ þ h1 q2 ; dt

ð1  b1 ÞG1  b2 G2  De3 þ De1 ˜ dD ˜ g þ h1 q3 ; þ h1 Dx dt dN˜ þ h1 q4 ; PP þ mDe2 þ De3 ¼ dt ¼

ð18Þ

at each time step. Here the qi are computed from Eq. (9). This set of equations imposes four constraints on six unknowns PP, G1, G2, Dei, i = 1,2,3 and we have to regularize it to get a unique solution. A natural

17

regularization is to seek the solution to Eq. (18) that would be closest to some first guess. Towards this end, we take parameterizations for PP, G1, G2, Dei, i = 1,2,3 and the optimized values of the model parameters and compute the first guess values P˜P, ˜ 1, G ˜ 2, D ˜ eii = 1,2,3 from the data. Then we solve the G minimization problem ðPP; G1 ; G2 ; G3 ; De1 ; De2 ; De3 Þ ˜ 1 Þ2 ˜ 2 þ ðG1  G ¼ Arg min½ðPP  PPÞ ˜ 2 Þ2 þ þ ðG2  G

3 X ˜ i Þ2  ðDei  De

ð19Þ

i¼1

subject to Eq. (18) at each time increment. Formally, we impose the equations (Eq. (18)) expressing the balance of mass between different biological compartments as strong constraints while the particular mathematical expressions used to describe the rates of corresponding processes are imposed as weak constraints. As a result, we obtain estimates of the fluxes (or, the rates of biological processes) that are compatible with the system evolution as provided by the data, at the cost of breaking the formulas used in the model to parameterize these processes. We note the difference with the strong constraint data assimilation where the estimates of the fluxes are computed from the model parameterizations, at the expense of getting a solution that does not fit the data within the data errors due to uncertainties in the parameterizations. After solving the secondary inversion for the fluxes, we can look at values of the terms in the right hand side of Eq. (19) and see which estimate deviates most notably from the first guess. This allows us to conclude which parameterization is mainly responsible for deviations of the model trajectory from the data, and thus which should be improved. The only point is that we have no measurements of all model state variables at each time increment. However, we can use the inverse solution obtained with the weak constraint data assimilation. This makes it possible to reduce the data misfits up to the level of the data errors. Fig. 8 depicts the estimates of the rates of biological processes obtained with the procedure described above for Bermuda Station ‘‘S’’. It is seen that PP and Gi, i = 1,2 are almost identical to P˜P and ˜ i, i = 1,2 and we can conclude that the primary G

18

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

Fig. 8. The seasonal evolution of the biological sources and sinks at Bermuda Station ‘‘S’’ as derived from the model solution (dashed line), the inverse solution (thin solid line) and from the secondary inversion (thick solid line). Stars denote the mean BATS primary production data.

production and the grazing rates are described in the model with good accuracy. The close fit of the curves to the measured primary production supports this conclusion. On the contrary, large deviations of Dei ˜ ei, i = 1,2,3 point out that the parameterizations from D of phytoplankton mortality De1, zooplankton losses De2 and detritus breakdown De3 contribute much to the model errors.

8. Conclusions When dealing with ecosystem models, we have to keep in mind that their structure is not completely derived from justified laws of the nature and thus the mathematical description of the biological processes is not universal. Thus, we are uncertain not only about values of the numerous model parameters, but also

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

about the model parameterizations and the model errors. In this situation, the weak constraint formulation provides an effective tool for tuning and calibrating ecosystem models. In this study, we have demonstrated two new outputs of weak constraint variational data assimilation which have been not considered previously and which are completely unrecoverable by a strong constraint method. The weak formulation of the model can shed light on the accuracy of the model parameterizations. We presented a procedure of the secondary inversion to estimate the biological fluxes, after getting the weak constraint inverse solution for the model compartments. For Bermuda Station ‘‘S’’, we found that the phytoplankton mortality, zooplankton losses and breakdown of detritus deviates notably from the first guess values, computed by means of the mathematical formulas which parameterize these processes. Consequently, these parameterizations make the major contribution to the model equation residuals obtained when minimizing the cost function (Eq. (3)) and thus should be improved. The procedure presented for the secondary inversion of fluxes makes it possible to restore the mass balance broken while performing the weak constraint parameter estimation and to refine the estimates of the biogeochemical fluxes. In addition, the data assimilation approach used in combination with a posteriori choosing of the data weight allowed us to compare the quality of the data and the model prediction. In particular, we were able to discern the low quality of the satellite data for high latitudes and for the coastal regions. The significant spatial variations of the optimized model parameters appeared to be in qualitative agreement with other estimates, implying that, generally, optimal values of model parameters obtained in a location cannot be successfully used for all other locations of the basin under consideration. Satellite ocean-color data are an extremely valuable source of information for studying the spatial variability of the model parameters. However, they may not suffice for discriminating between different system trajectories resulting in similar chlorophyll concentrations, especially if the data are not accurate enough and a ‘‘combination with in situ observations is vital’’ (Hemmings et al., 2003).

19

These issues are not new and only confirm some notions existing in the ocean biogeochemical community. However, the fact of this confirmation per se means that the data assimilation approach used have passed a good test and can be applied to more sophisticated models. Returning back to the two questions – are these estimates applicable at other locations? – are the model errors truly negligible? formulated in Section 1, we should say that our answers are negative, and thus, the problem of tuning ecosystem models seems to become more challenging than it has been considered before.

Acknowledgements The authors are grateful to the editors and anonymous reviewers whose comments were very helpful and allowed us to improve the paper and to our colleague Markus Schartau (AWI), discussions with whom were of great value.

References Antoine, D., Andre, J.-M., Morel, A., 1996. Oceanic primary production: II. Estimation at global scale from satellite (Coastal Zone Color Scanner) chlorophyll. Glob. Biogeochem. Cycles 10, 57 – 69. Balch, W., Evans, R., Brown, J., Feldman, G., McClain, C., Esaias, W., 1992. The remote sensing of Ocean Primary Productivity: use of a new data compilation to test satellite algorithms. J. Geophys. Res. 97 (C2), 2279 – 2293. Bennett, A.F., 1992. Inverse Methods in Physical Oceanography. Cambridge Univ. Press, Cambridge, UK. 346 pp. Bennett, A.F., Thorburn, M.A., 1992. The generalized inverse of a nonlinear quasigeostrophic model. J. Phys. Oceanogr. 22, 213 – 230. Brummelhuis, P.G.J., Heemink, A.W., van den Boogaard, H.F.P., 1993. Identification of shallow sea models. Int. J. Numer. Methods Fluids 17, 637 – 665. Cloern, J.E., Grenz, C., Vidergar-Lucas, L., 1995. An empirical model of the phytoplankton chlorophyll: carbon ratio—the conversion factor between productivity and growth rate. Limnol. Oceanogr. 40, 1313 – 1321. Conkright, M.E., Levitus, S., Boyer, T.P., 1994. World Ocean Atlas. Nutrients, vol. 1. NOAA, Washington, DC. 150 pp. Cramer, H., 1954. Mathematical Methods of Statistics, 6th ed. Princeton Univ. Press, Princeton, N.J., USA. Eknes, M., Evensen, G., 1997. Parameter estimation solving a weak

20

S.N. Losa et al. / Journal of Marine Systems 45 (2004) 1–20

constraint variational formulation for an Ekman model. J. Geophys. Res. 102, 12479 – 12491. Eknes, M., Evensen, G., 2002. An ensemble Kalman filter with a 1-D marine ecosystem model. J. Mar. Syst. 36, 75 – 100. Evans, G.T., 1999. The role of local models and data sets in the Joint Global Ocean Flux Study. Deep-Sea Res., Part 1, Oceanogr. Res. Pap. 46, 1369 – 1389. Evensen, G., Fabrio, N., 1997. Solving for the generalized inverse of the Lorenz model. J. Meteorol. Soc. Jpn. 75, 229 – 243. Evans, G.T., Parslow, J.S., 1985. A model of annual plankton cycles. Biol. Oceanogr. 3, 327 – 347. Evans, R.H., Gordon, H.R., 1994. Coastal Zone Color Scanner ‘‘system calibration’’: a retrospective examination. J. Geophys. Res. 99, 7293 – 7307. Fasham, M.J.R., 1993. Modelling the marine biota. In: Heimann, M. (Ed.), The Global Carbon Cycle. Springer-Verlag, New York, pp. 504 – 547. Fasham, M.J.R., Evans, G.T., 1995. The use of optimization techniques to model marine ecosystem dynamics at the JGOFS station at 47jN 20jW. Philos. Trans. R. Soc. Lond., B 348 (1324), 203 – 210. Fennel, K., Losch, M., Schroeter, J., Wenzel, M., 2001. Testing a marine ecosystem model: sensitivity analysis and parameter optimization. J. Mar. Syst. 28, 45 – 63. Geider, R.J., Macintyre, H.L., Kana, T.M., 1996. A dynamic model of photoadaptation in phytoplankton. Limnol. Oceanogr. 41, 1 – 15. Ghil, M., Malanotte-Rizzoli, P., 1991. Data assimilation in meteorology and oceanography. Adv. Geophys. 33, 141 – 266. Gong, J., Wahba, J., Johnson, D.R., Tribbia, J., 1998. Adaptive tuning of numerical weather prediction models: simultaneous estimation of weighting, smoothing and physical parameters. Mon. Weather Rev. 125, 210 – 231. Hemmings, J.C.P., Srokosz, M.A., Challenor, P., Fasham, M.J.R., 2003. Assimilating satellite ocean-colour observations into oceanic ecosystem models. Philos. Trans. R. Soc. Lond. A361, 33 – 39. Hurtt, G.C., Amstromg, R.A., 1996. A pelagic ecosystem model calibrated with BATS data. Deep-Sea Res. 43, 653 – 683. JGOFS Rep., 1997. In: Evans, G.T., Garcon, V. (Eds.), One dimensional models of water column biogeochemistry. JGOFS Report No. 23. Feb. 1997. p. 85. Kagan, B.A., Kivman, G.A., Ryabchenko, V.A., Losa, S.N., 1997. Assimilation of satellite chlorophyll ‘a’ concentrations into an ecosystem model of the ocean upper mixed layer. Dokl. Ross. Akad. Nauk. 355, 688 – 692 (in Russian). Kivman, G.A., Kurapov, A.L., Guessen, A.V., 2001. An entropy approach to tuning weights and smoothing in the generalized inversion. J. Atmos. Ocean. Technol. 18, 266 – 276. Kyewalyanga, M.N., Platt, T., Sathyendranath, S., Lutz, V.A., Stuart, V., 1998. Seasonal variations in physiological parameters of phytoplankton across the North Atlantic. J. Plankton Res. 20, 17 – 42. Levitus, S., Boyer, T.P., 1994. World Ocean Atlas. Temperature, vol. 4. NOAA, Washington, DC. 117 pp.

Matear, R.J., 1995. Parameter optimization and analysis of ecosystem models using simulated annealing: a case study at Station P. J. Mar. Res. 53, 571 – 607. Natvik, L.-J., Eknes, M., Evensen, G., 2001. A weak constraint inverse for a zero dimensional ecosystem model. J. Mar. Syst. 28, 19 – 44. Oberhuber, J.M., 1988. An Atlas Based on the COADS Data Set: The Budgets of Heat, Buoyancy and Turbulent Kinetic Energy at the Surface of the Global Ocean. Max Plank Institute for Meteorology, Hamburg, Germany. Platt, T., Caverhill, C., Sathyendranath, S., 1991. Basin-scale estimates of oceanic primary production by remote sensing: the North Atlantic. J. Geophys. Res. 96, 15147 – 15159. Popova, E.E., 1995. Non-universal sensitivity of a robust ecosystem model of the ocean upper mixed layer. Ocean Model. 109, 2 – 5. Popova, E.E., Fasham, M.J.R., Osipov, A.V., Ryabchenko, V.A., 1997. Chaotic behaviour of an ocean ecosystem model under seasonal external forcing. J. Plankton Res. 19, 1495 – 1515. Prunet, P., Minster, J.-F., Ruiz-Pino, D., 1996a. Assimilation of surface data in a one-dimensional physical – biogeochemical model of the surface ocean: 1. Method and preliminary results. Glob. Biogeochem. Cycles 10, 111 – 138. Prunet, P., Minster, J.-F., Ruiz-Pino, D., 1996b. Assimilation of surface data in a one-dimensional physical – biogeochemical model of the surface ocean: 2. Adjusting a simple trophic model to chlorophyll, temperature, nitrate, and pCO2 data. Global Biogeochem. Cycles 10, 139 – 158. Schartau, M., Oschlies, A., Willebrand, J., 2001. Parameter estimates of a zero-dimensional ecosystem model applying the adjoint method. Deep-Sea Res., Part 2, Top. Stud. Oceanogr. 48 (2001), 1769 – 1800. Srokosz, M.A., 1998. Data assimilation into oceanic ecosystem models: a critical review, some comparisons, and recommendations. Southampt. Oceanogr. Cent. Intern. Doc. 32 32 pp. (available at: http://www.soc.soton.ac.uk/JRD/SAT/BlueSkies/SOC_ Report_32/SOC_No_32.pdf). Spitz, Y., Moisan, J.R., Abbott, M.R., Richman, J.G., 1998. Data assimilation and a pelagic ecosystem model: parameterization using time series observations. J. Mar. Syst. 16, 51 – 68. Spitz, Y., Moisan, J.R., Abbott, M.R., 2001. Configuring an ecosystem model using data from the Bermuda Atlantic Time Series (BATS). Deep-Sea Res. 48, 1733 – 1768. Vinogradov, M.E., Shushkina, E.A., Vedernikov, V.I., Gagarin, V.I., Nezlin, N.P., Shebertsov, S.V., 1995. Characteristics of the Pacific epipelagic ecosystems based on the satellite and expedition data: abiotic parameters, and production indices of phytoplankton. Oceanolgy, English Translation 35, 208 – 217. Tarantola, A., 1987. Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation Elsevier, Amsterdam. 613 pp. Wahba, G., 1990. Spline Methods for Observational Data. SIAM, Philadelphia. 169 pp. World Ocean Database, 2001. http://www.nodc.noaa.gov/OC5/ WOD01/pr_wod01.html.