Simultaneous data-based optimization of a 1D-ecosystem model at three locations in the North Atlantic:

Simultaneous data-based optimization of a 1D-ecosystem model at three locations in the North Atlantic:

Part I—Method and parameter estimates

by Markus Schartau^1,2and Andreas Oschlies³ ABSTRACT

An optimization experiment is performed with a vertically resolved, nitrogen-based ecosystem model, composed of four state variables (NPZD-model): dissolved inorganic nitrogen (N), phyto- plankton (P), herbivorous zooplankton (Z) and detritus (D). Parameter values of theNPZD-model are optimized while assimilating observationsat three locations in the North Atlantic simultaneously, namely at the sites of the Bermuda Atlantic Time-Series Study (BATS; 31N 64W), of the North Atlantic Bloom Experiment (NABE; 47N 20W), and of Ocean Weather Ship-India (OWS-INDIA;

59N 19W). A method is described for a simultaneous optimization which effectively merges different types of observational data at distinct sites in the ocean. A micro-genetic algorithm is applied for the minimization of a weighted least square misŽ t function. The optimal parameter estimates are shown to represent a compromise among local parameter estimates that would be obtained from single-site optimizations at the individual locations. The optimization yields a high estimate of the initial slope parameter of photosynthesis(a), which is shown to be necessaryto match the initial phases of phytoplankton growth. The estimate ofa is well constrained by chlorophyll observations at the BATS and OWS-INDIA sites and likely compensates for a deŽ ciency in the parameterizationof light-limited growth. The optimization also points toward an enhanced recycling of organic nitrogen which is perceived from a high estimate for the phytoplankton mortality/

excretion rate.

1. Introduction

One general task of marine ecosystem models is to represent those ecological processes which signiŽ cantly contribute to biogeochemical  uxes in the ocean. For this purpose, a variety of biological models have been developed (e.g. Evans and Parslow, 1985; Fasham et al.,1990; Steele and Henderson, 1992; Hurtt and Armstrong, 1999; Doneyet al.,1996;

Mooreet al.,2001). These models differ in complexity, from simple models containing three biological state variables up to more complex ones with, presently, some thirty compartments. While more complex models appear to be more realistic at Ž rst glance, they

1. Alfred-Wegener-Institut fu¨r Polar- und Meeresforschung, Bussestr. 24, 27670 Bremerhaven, Germany.

2. Present address: Marine Sciences Research Center, SUNY, Stony Brook, New York, 11794, U.S.A.email:

markus.schartau@stonybrook.edu

3. Institut fu¨r Meereskunde, Du¨sternbrooker Weg 20, 24105 Kiel, Germany.

765

usually contain a large number of parameterizations with poorly known parameters that are often difŽ cult to be determined directly from observational or physiological information.

As soon as a model becomes too complex, in the sense that it contains unconstrainable parameters (degrees of freedom), it may not be useful for extrapolation to climatic conditions other than those used to tune the model. Hence, the available data places signiŽ cant limitations on the necessary model complexity (Matear, 1995). When, on the other hand, the model is too simplistic, the simulated biogeochemical  uxes might be seriously wrong due to the absence of important processes. Yet, the minimum complexity that is required for an ecosystem model to reliably reproduce biogeochemical observations while being fully constrainable has not been explored. We believe that an objective data assimilation method will help to resolve this issue. The present study provides tools that may subsequently be applied to a more systematic assessment of different marine ecosystem models, a task far beyond the scope of this work.

The nitrogen-based ecosystem model developed by Fasham et al. (1990) (hereafter named FDM-model) has become a standard model used in various studies ranging from zero-dimensional mixed-layer applications to fully three-dimensional coupled ecosystem- circulation models. Regarding basin-scale coupled physical-biological simulations, as presented by Sarmientoet al.(1993) and Fashamet al.(1993), the FDM-model is one of the more complex models implemented so far. Some of its parameter values were adopted from published laboratory experiments whereas others were approximated from rates derived fromin situmeasurements. During the subsequent years, the need to seek optimal parameter estimates for the FDM-model from data assimilative investigations was stressed repeatedly (Fashamet al.,1993; Fasham and Evans, 1995).

Local parameter optimizations of the FDM-model were performed at the US-JGOFS station of the Bermuda Atlantic Time-Series Study (BATS) near 31N, 64W and at 47N, 20W (site of the North Atlantic Bloom Experiment, NABE) by Spitzet al. (1998) and Fasham and Evans (1995), respectively. In both studies a zero-dimensional model of the upper ocean’s mixed layer was used. Fasham and Evans (1995) could reproduce the observed nitrate and chlorophyll concentrations at the NABE site, as well as the primary production with an optimized set of parameter values. At the BATS site, however, data-assimilation experiments with the FDM-model seemed to be more problematic (Spitz et al., 1998, 2001). These studies also showed that, despite the large number of high- quality data available near Bermuda, not all parameters could be constrained by the observations, suggesting that even a 7-component model may have too many degrees of freedom with respect to Ž tting well-sampled time series data.

Developing an elaborate and simpliŽ ed FDM-model with a reduced number of parame- ters, all of which could be well constrained, was the motivation of Hurtt and Armstrong (1996). They proposed a zero-dimensional model which accounts for allometric relation- ships for the biological rate parameters. A modiŽ ed version (hereafter named HA-model) was then optimized to data at two locations in the North Atlantic simultaneously (Hurtt and Armstrong, 1999). To simultaneously Ž t observations at the two sites, namely BATS and

the location of the Ocean Weather Ship India (OWS-INDIA), the authors had to include a parameterization of iron limitation in their model, which they assumed to be effective at OWS-INDIA but not at the BATS location. To our knowledge, this promising model approach has not yet been applied to a basin-wide 3D-model. A different series of model studies was performed with another reduced FDM-model in a coupled physical-biological model of the North Atlantic (Oschlies and Garc¸on, 1999; Oschlies et al., 2000). Their ecosystem model combined nitrate and ammonium to dissolved inorganic nitrogen (N).

Model compartments such as phytoplankton (P), zooplankton (Z) and detritus (D) remained in the model (hence referred to asNPZD-model) whereas bacteria and dissolved organic nitrogen were not resolved explicitly. For the coupled simulations they relied on parameter values similar to those published by Sarmientoet al.(1993) and Fashamet al.

(1993).

The main objective of this study is to identify a single set of parameter values that improves the performance of the NPZD-model at three different locations in the North Atlantic where time-series data are available: at the BATS site (31N, 64W), at the NABE site (47N, 20W), and at OWS-INDIA (59N, 19W) (Fig. 1). We attempt to provide optimal parameter estimates for theNPZD-model which can subsequently be used in a basin-scale simulation of the North Atlantic. This is achieved by assimilating observations which were Figure 1. North Atlantic Ocean. The three locations of BATS (31N 64W), of NABE (47N 20W) and of OWS-INDIA (59N 19W) used in our data-assimilative investigationsare marked. Gray shaded contours show ocean bathymetry.

collected mainly as part of the Joint Global Ocean Flux Study (JGOFS) at the three sites. In the present paper we investigate optimal parameter estimates and the applicability and robustness of the optimization procedure. Particularly, we will discuss which parameter estimates are likely to compensate for deŽ ciencies of theNPZD-model. A detailed analysis of standing stocks and biological  uxes simulated by the optimized model is the subject of an accompanying paper (Schartau and Oschlies, 2003; this issue).

2. Method

a. Model description

The model’s biological state variables (C_i) are strongly simpliŽ ed representations of either nutrients or organisms that are assumed to be evenly spread within a grid box. The vertical distribution of the state variables is simulated as a function of time, with turbulent vertical mixing coefŽ cients taken from a 3D physical model run (see below). Effects due to vertical advective  ow, and horizontal  uxes resulting from a divergent vertical  ow Ž eld, are not accounted for. The governing one-dimensional equation can be formulated as follows:

]Ci

]t 52wi

]Ci

]z 1 ]

]z

S

D

K_r is the turbulent mixing coefŽ cient and w_i is the sinking velocity which becomes nonzero only for detritus. The terms on the right-hand side represent sinking of detritus, turbulent mixing and source minus sink (SMS) terms that describe the inherent biological processes.

The biological interactions among the four compartments of theNPZD-model are sketched in Figure 2. The arrows show the nitrogen  uxes, with symbols indicating those parameters that are associated with the rates for each particular  ux. In contrast to previous model versions of thisNPZD-model (Oschlies and Garc¸on, 1999; Oschlies, 2001), not only phytoplankton growth but also all remineralization rates, i.e. the  uxes from P,Z, andDtoN, are temperature dependent. All model parameters are listed in Table 1. The full model equations are listed in the Appendix. The biological model is initialized with vertical nitrate proŽ les from Conkrightet al.(1994). The time-steps of integration are 15 minutes for the biological state variables whereas they are 1 hour for the advective-diffusive equation. The spin-up time equals two identical years with the same physical components as derived for the year 1989.

b. Physical forcing

With the aim to achieve a high degree of consistency with a basin-scale ecosystem- circulation model, the one-dimensional ecosystem model is embedded into a physical environment taken from three-dimensional ocean circulation simulations. The chosen

ocean circulation model is identical to the one applied by Oschlies and Garc¸on (1999). The horizontal resolution allows for mesoscale variability (“eddy permitting”) with a meridi- onal grid of 1/3 and a zonal grid spacing of 2/5. The water column is partitioned into 37 levels with the Ž rst 10 levels resolving the upper 126 meter (see Table 2). Resolving the entire water column down to the bottom excludes the need for an open lower boundary of the model. The fact that vertical grid spacing increases from about 30 m below the euphotic zone to 250 m below 1000 m is considered to have only minor impact on the simulated nitrate supply by deep winter mixing, because nutrient gradients at depth are usually very weak. On the few-year time scales considered here, interactions with the sediment can be neglected.

The 3D-model run from which the physical environment is taken at the three 1D sites was forced with daily mean reanalysis data of the European Center of Medium Range Weather Forecast (ECMWF), covering the period 1989 through 1993 (Nicolas Ferry, personal communication). The one-dimensional ecosystem simulations run here use the same ECMWF surface short-wave radiation data as the basin-scale simulation.

Figure 2. Structure of the ecosystem model. The compartments (state variables) are dissolved inorganic nitrogen (N), phytoplankton biomass (P), herbivorous zooplankton (Z) and detritus (D).

The arrows indicate the direction of mass  ux. Those parameters which control a particular mass

 ux are additionally listed. The parameter symbols are explained in Table 1.

c. Observations

At all three locations Ž ve types of observations are considered in this study: nitrate1 nitrite (DIN), chlorophylla (Chla), ¹⁴C-primary production (C-PP), particulate organic nitrogen (PON) and zooplankton biomass (ZOO). Whenever nitrite measurements are not available, then only nitrate is considered. All observations are Ž rst interpolated onto a 1 meter vertical grid and then averaged over the model’s grid boxes. The average of those individual proŽ les that are available for each month are calculated for each data type within the upper Ž fteen model layers, covering a depth range of 411 meters which is usually well resolved by the measurements. Note that these monthly-mean values have to be derived from a relatively sparse sampling schedule (particularly at the NABE and OWS-INDIA sites with sometimes none or only a single proŽ le per month). Hence, these monthly values do not necessarily represent true monthly means. Table 3 gives an overview of the observations considered.

Location 31N 64W.Near Bermuda, data are available from the BATS, as a part of the U.S.

JGOFS project (Michaels and Knap, 1996). The BATS data are provided by the Bermuda Biological Station for Research (BBSR). [BATS extraction sitehttp://www.bbsr.edu/users/

ctd/] Except for zooplankton biomass, all monthly mean values are obtained from Table 1. Parameters of the NPZD-model. The Ž rst thirteen parameters of the list enter the

optimization process.

p_n Symbol Parameters for variation Unit

p₁ m_m Growth rate parameter d²¹

p₂ a Slope of photosynthesisversus light intensity

m²W²¹d²¹

p₃ F_m^p Phytoplankton loss rate parameter d²¹

p₄ F*_p Phytoplankton quadratic loss m³mmol N²¹d²¹

p₅ k_N Half saturation constant ofNuptake rate mmol Nm²³ p₆ k Attenuation coefŽ cient due to phytoplankton m²mmol N²¹

p₇ g Maximum grazing rate d²¹

p₈ e Prey capture rate m⁶mmol N²²d²¹

p₉ F_m^z Herbivore loss rate parameter d²¹

p_{1 0} F*_z Herbivore quadratic mortality m³mmol N²¹d²¹ p_{1 1} b Assimilation efŽ ciency of herbivores dimensionless p_{1 2} g_m Remineralization rate parameter of detritus d²¹

p_{1 3} w_s Detrital sinking velocity m d²¹

Fixed parameters Value and unit

k_w Light attenuation due to water 0.04 m²¹

R molar carbon to nitrogen ratio 6.625

f_{P A R} short-wave fraction of PAR 0.43

C_{r e f} growth coefŽ cient 1.066

c growth coefŽ cient 1.000 (°C)²¹

Table 2. Vertical levels of the numerical model grid. Units are meters.

Vertical layers of theNPZD-model

Model level Depth of grid point

Depth of grid

box bottom Thickness of grid box

1 5.50 11.00 11.00

2 17.00 23.00 12.00

3 29.00 35.00 12.00

4 41.00 47.00 12.00

5 53.00 59.00 12.00

6 65.50 72.00 13.00

7 78.50 85.00 13.00

8 91.50 98.00 13.00

9 104.50 111.00 13.00

10 118.50 126.00 15.00

11 140.50 155.00 29.00

12 179.55 204.09 49.09

13 232.60 261.10 57.01

14 295.03 328.95 67.85

15 370.21 411.47 82.52

16 462.51 513.54 102.07

17 577.37 641.19 127.65

18 721.47 801.74 160.55

19 900.89 1000.04 198.30

20 1125.04 1250.04 250.00

21 1375.04 1500.04 250.00

22 1625.04 1750.04 250.00

23 1875.04 2000.04 250.00

24 2125.04 2250.04 250.00

25 2375.04 2500.04 250.00

26 2625.04 2750.04 250.00

27 2875.04 3000.04 250.00

28 3125.04 3250.04 250.00

29 3375.04 3500.04 250.00

30 3625.04 3750.04 250.00

31 3875.04 4000.04 250.00

32 4125.04 4250.04 250.00

33 4375.04 4500.04 250.00

34 4625.04 4750.04 250.00

35 4875.04 5000.04 250.00

36 5125.04 5250.04 250.00

37 5375.04 5500.04 250.00

bi-weekly to monthly data, covering the years 1989–1993. Zooplankton biomass is taken from measurements of heterotrophic nano- and microzooplankton (Caronet al.,1995).

47N 20W.Here, most measurements were taken during NABE (Ducklow and Harris, 1993) in the year 1989. Therefore, the monthly observational representatives are biased toward 1989. The initial data base was received from the British Oceanographic Data Center (BODC) and extended by German JGOFS investigations until 1996 (Koeve, personal communication), also available through the German JGOFS data management. All data entering the calculations are selected from a 53 5 degree area (17.5W–22.5W, 44.5N–

49.5N). Microzooplankton biomass data are taken from Fasham and Evans (2000) who referred to measurements of Verityet al.(1993) and two additional observations made in late summer of 1989.

Table 3. Collected data for the data-assimilationexperiment at three locations in the North Atlantic.

Data types: dissolved inorganic nitrogen (DIN); Chlorophylla (Chl); Carbon-based Primary Production (CPP); Particulate Organic Nitrogen (PON); Zooplankton biomass given in nitrogen units (ZOO). Data sources:¹US-JGOFS, BBSR;²Caronet al.(1995);³Koeve, personal communi- cation;⁴Roy Lowry (BODC), personal communication;⁵Fasham and Evans (2000), Verityet al.

(1993), and Fasham, personal communication;⁶Williams (1988), and data from 1996–1997 collected as part of PRIME. The number of proŽ les used for the construction of monthly representatives,and the seasonal coverage of the Ž nal climatology.

Area of observation Data type

Number of

proŽ les Years Seasonal representation/month/

31N 20W (BATS site)

DIN¹ 77 1989–1993 12 months

Chl¹ 75 1989–1993 12 months

CPP¹ 67 1989–1993 12 months

PON¹ 77 1989–1993 12 months

ZOO² 6 1989, 1990 3 months; 3/4/8

47N 20W (NABE site)

DIN³ 96 1989–1997 9 months; 1/2/3/4/5/7/8/9/10 Chl³ 100 1989–1997 9 months; 1/2/3/4/5/7/8/9/10

CPP⁴ 29 1989 4 months; 4/5/7/8

PON⁴ 48 1989–1993 3 months; 4/5/7

ZOO⁵ 12 1989 3 months; 5/6/7

59N 19W

(OWS-INDIA site)

DIN⁶ 97 1971–1975,

1996

8 months; 3/4/5/6/7/8/9/10

Chl⁶ 335 1971–1975,

1996–1997

8 months; 3/4/5/6/7/8/9/10

CPP⁶ 78 1970–1975,

1997

8 months; 3/4/5/6/7/8/9/10

PON⁶ 14 1996 2 months; 6/7

ZOO⁶ 20 1996 2 months; 6/7

59N 19W.Most observations yielding a good seasonal coverage at this northern Atlantic site were taken during the years 1971–1974 (Williams, 1988). Although some additional data were collected in 1989, 1996, and 1997 as part of the Plankton Reactivity In the Marine Environment (PRIME) project; the derived monthly means are dominated by data from the early 1970s. As for the NABE site, only data within a 5 3 5 degree area are considered (16.5W–21.5W, 56.5N– 61.5N).

d. Cost function

i. DeŽ nition. The cost function is deŽ ned as a sum of weighted least square misŽ ts between model results and observations. In order to compute these misŽ ts one could sample the model at the same times and locations as the observations were taken. However, with the physics taken from an eddy-permitting model forced by daily atmospheric winds and heat  uxes, large model-data misŽ ts can arise already from small phase errors in the model. For example, misplaced eddies or systematic errors in the buoyancy budget may easily change the onset of the spring bloom by several days to weeks. Temporal weighting terms could be introduced but are hardly justiŽ able when assimilating observations at OWS-INDIA, which were predominantly taken in periods (1971–74) not covered by the model forcing (1989–93). To reduce the likely mapping of physical phase errors into the cost function and hence into the optimization of biological model parameters, it was decided to use only monthly-averaged data and monthly-averaged model results instead.

Because the available observations do not resolve interannual variability of the simulation period at OWS-INDIA and only marginally at the NABE site, it was further decided to map all available observations of all years into a single “climatological” composit for each month. Initial tests using a cost function composed of misŽ ts of such monthly “climatological”

observations and “climatological” monthly model means were not successful because the cost function turned out to be insensitive to unrealistic model drifts over the Ž ve-year simulation period. In order to get a better handle on such model drifts, we proceeded by comparing the simulated monthly means of each individual year of the Ž ve-year simulation with the monthly

“climatological” observations. Still, results were not fully satisfactory, and eventually we resorted to include an additional steady-state constraint to the cost function (see below).

At the end of this iterative process, the total cost function used for the simultaneous optimization at the three sites is deŽ ned as follows: the total cost combines the individual cost function contributions,7_l, from the three locations:

7total51 2

O

l51 3

7l (2)

withlbeing the location index. The individual cost functions are split up into two parts:

7_l5

O

y51 5 7yl

Cl

1 1

s²~N₈₉^total2N_y^total!l2 (3)

with

Nytotal5

E

E

The Ž rst term describes the data-model misŽ ts and includes a scaling factor C_lfor the individual stations discussed below. The second term is the steady-state constraint which penalizes large deviations from the total nitrogen inventory reached after spin up. The integralN_y^totalstands for the total nitrogen mass within the upper 411 meters of the yeary.

This integral is calculated for each year and compared with the initial nitrogen inventory (after the model’s two-year spin up) of the year 1989 (N₈₉^total). This second term is based on a plausible, albeit subjective argument, in order to avoid large systematic drifts that can occur in the simulated nitrogen inventory at the individual stations. It actually states that biological production, based on newly entrained nutrients (new production), should not deviate too much from the biogenic export  ux. Note, however, that this term will also penalize interannual variability in the model as far as the total nitrogen inventory is affected. The corresponding standard deviationsis chosen to be 5% of the initial inventory of 1989. At every location the annual data-model misŽ t contributionsto the cost function are calculated as follows:

7_yl5

O

j

J 1

1_jl

O

m 1jl

O

K 1

sj2~f_j2y_j^obs!mk2 (5) with the monthly “climatological” representativey_jm^obsof typejand the modeled monthly mean (f_jm) of yeary. The total number of observational types isJ55 (DIN, Chla, C-PP, PON, and ZOO).1_jlis the number of months for which observations of typejat locationl are available. The maximal depth for data to be assimilated into the model is 411 meters which is equivalent to the bottom of theK515 grid box. Correlations among the different variables are not accounted for. Division by the number of available observational months gives same weights to observations obtained only during a few months (e.g. zooplankton) compared to those covering a more complete seasonal cycle, see Table 3. The assigned weights for DIN ares_DIN50.1 mmol N m²³. For Chlaand PON concentrations weights ofs_Chl5 0.01 mg m²³ands_PON5 0.0357 mmol m²³(’0.5 mg m²³) are assumed, respectively. Because of large differences in primary production among the three sites, different weights of 15% of the annual mean primary production at the respective site are considered, yieldingsPP50.8 mg C m²³d²¹for the OWS-INDIA site, 1.0 mg C m²³d²¹ for the NABE site and 0.3 mg C m²³d²¹for the BATS site. The weights for zooplankton biomass are set tos_ZOO50.01 mmol N m²³at all locations.

ii. Scaling of stations.The simplest idea for an overall cost function is to add together all weighted least square misŽ ts of the three locations. Unfortunately, this results in a strong bias of the solution toward one particular location. For instance, as soon as some sort of observation (e.g. chlorophyll concentrations)differs between two locations by one order of

magnitude, the same relative model-data misŽ ts produce different absolute contributions to the cost function. In test experiments without any scaling of the locations, the optimization converged toward parameter values that generated a model solution with extinct biology at the location with lowest productivity and biomass (in our case at the BATS site). To ensure that similar relative misŽ ts at the different locations give similar cost function contribu- tions, a scaling factorC_lis introduced that considers averaged observational values at the respective locationl:

Cl51 2

O

j51 J ~y¯^obs!_j²

sj2 (6)

with the subscriptjreferring to the data type. At each location the square of the averaged observation (in space and time) is divided by the assumed weights, as prescribed above.

e. Optimization procedure

i. Micro-genetic algorithm (mGA) for optimization.Schartau et al. (2001) showed that under realistic conditions, especially if model deŽ ciencies exist, a gradient technique can produce optimal solutions which are sensitive to the initial parameter guesses. For robust parameter estimates they had to perform hundreds of individual optimizations, starting from a variety of initial parameter guesses. Intercomparisons of various methods for the optimization of parameter values of marine biogeochemical models (Vallino, 2000; Athias et al., 2000) revealed some advantages of sophisticated stochastic algorithms such as genetic algorithms. Stochastic or quasi-stochastic algorithms do not require special programming efforts for the calculation of the cost function’s gradient (e.g., coding of an adjoint model). Furthermore, they are more robust than pure gradient-descent techniques in cases when the cost function contains regions with plane-geometry that result from low sensitivities to parameter variations. For this reason, we resort to the concept of genetic algorithms (GA) (Holland, 1975; Goldberg, 1989).

When applying a GA for parameter optimization, a single set of parameters (parameter vector) is represented by an individual which is coded as a binary string (similar to genes of a chromosome). A generation is set up by a prescribed number (n) of individuals.Selection of individuals, recombination of genes (crossover), and mutation of individuals are the basic operations regarded in genetic algorithms. The selection process follows the evolutionary principle of “survival of the Ž ttest,” with the Ž tness being expressed by a small cost function value. The recombination operation describes the exchange of genes among the selected individuals (parents), setting up an offspring generation with new parameter vectors (children). Often mutation is also regarded. Mutation induces small modiŽ cations to some children and therefore brings in new information that is independent of the selected parents.

In this study, we apply a micro-genetic algorithm (mGA), coded and published by Carroll (1996). Numerical details of themGA are described by Krishnakumar (1989) who

Ž rst presented the algorithm as being able to optimize nonstationary functions. ThemGA is based on the same operations as the general GA, but it does not contain mutation and gives greater emphasis on elitism principles. An elitism operator assures that information of the very best individual (parameter set) is retained for recombination. As soon as all individuals of one generation show less than 5% difference among each other, then a new random population is generated, apart from the best individual which is saved. This process is repeated several times. Hence, while general convergence is achieved, the full parameter space is explored further. This convergence characteristic is well suited for our optimiza- tion problem.

The population size is set equal to the numbers of parameters of interest,n513. This is not mandatory, but was found to be a good choice in test experiments (Schartau, 2001). A redundant control parameter is added to the optimization process. Such a control parameter has no effect on the cost function and thereby allows testing for erroneous convergence.

The total number of generations is set to 2000, for which all experiments showed a well-converged cost function value. For recombination a single-pointed crossover (ex- change of bit-strings) is applied with a crossover probability of 1.0. Table 4 lists the different components of the parameter vector to be optimized, their respective upper and lower bonds and the resolution used by the optimization algorithm to generate random samples in parameter space.

ii. Errors of parameter estimates. In order to assess the reliability of the optimization algorithm and to obtain an estimate of the errors of the respective optimal parameters Table 4. Parameter setup for optimization. One set of parameter values (parameter vector) is representsby an individualwith a binary string. The increments yield the highest precision that, for each parameter, can be achieved with themGA. The length of a single bit string of one parameter (or number of binary digits) describes the number of possible values between lower and upper bound.

Parameter conŽ guration for micro-genetic algorithm

Parameter Lower bound Upper bound Increment # of possibilities

m_m 0.200 1.470 0.010 128

a 0.001 0.256 0.001 256

F_m^p 0.000 0.635 0.005 128

F*_p 0.010 0.955 0.015 64

k_N 0.100 0.730 0.010 64

k 0.010 0.073 0.001 64

g 0.025 1.600 0.025 64

e 0.025 1.600 0.025 64

F_m^z 0.000 0.635 0.005 128

F*_z 0.010 0.955 0.015 64

b 0.300 0.935 0.005 128

g_m 0.020 0.146 0.002 64

w_s 1.000 128.0 1.000 128

generated by this method, a total of four optimizations are performed (Fig. 3). The basic concept is to run additional optimizations with the original data being replaced by a synthetic data set,y^syn, generated by adding Gaussian noisehto the observations, with a noise variance equivalent to the error variancess_j²as used in the cost function Eq. 5:

yjkmsyn;yjkmobs1h~sjobs!, (7) withj,k,mreferring to data type, depth level, and month, respectively. This approach is related to the Monte Carlo bootstrap method (Efron, 1994). Computational costs restrict us to only three synthetic data sets that are utilized for additional optimizations,of which each single optimization requires 26000 model runs per location. On a 677 MHz DEC/ALPHA work station this results in a total of 308 hours CPU time for the entire procedure.

The optimization with the original data sety_okm^obs produces a best estimate,pˆ₀, of the parameter vector. Together with the three parameter vectors pˆ_i^syn that result from the optimizations with the three synthetic data setsy_jkn^syn, four independent optimization results (realizations) are then available. An estimate for the expected error in the optimal parameter setpˆ₀can then be constructed from two terms: A Ž rst error estimate,~sˆsd!n, is given by the standard deviation with respect to the mean value obtained from all four

Figure 3. Sketch of the optimization procedure.

realizations. A second term, (eˆ)_naccounts for the deviation of the best parameter estimate pˆ₀from the mean value^pˆ&:

~eˆ!n5~pˆ₀2^pˆ&!n. (8) Assuming these two error contributions to be independent, the approximate error variance of the optimal parameter set becomes

~sˆ!n5

Î

wheren refers to the respective component of the parameter vector. Note that this error estimate includes parameter uncertainties that arise from errors in the individual observa- tions, from possible redundancies in the model structure that give rise to parameter combinations that cannot be constrained by any data, and from potential deŽ ciencies of the optimization routine that may induce some scatter in the Ž nal solutions.

3. Results

a. Minimization of the cost function

All four optimizations performed (real observations plus three sets of observations with added noise) result in a signiŽ cant decrease of the cost function values (hereafter called

‘costs’) with respect to the initial values. The initial costs of the Ž rst generation range between7_total5157 and 1503, whereas the minimum yields a cost of7_total561 (a 61%

decrease for the lower end of the range of initial values). Since no information on the probability density function exists with regard to our optimization problem, it is impossible to describe an expectation value for the cost function of an optimized model that is consistent with the observations. Nevertheless, a reference cost function value can be derived by generating pseudo observations from the optimized model. These model data are sampled and processed in the same way as the real observations. By construction, these pseudo data are fully consistent with the model. When computing the cost function for these pseudo data in the same way as described for the real data in the method section, the resulting reference cost accounts for all preprocessing errors, e.g. those resulting from biases in the monthly data averages due to sparse sampling. Even a perfect model could not produce a cost function value lower than that given by this reference value. In Table 5 the cost function contributions at all three locations are shown for the optimal conŽ gurations together with the reference values. In addition, costs of a traditional parameter conŽ gura- tion are presented which rely on parameter values similar to those applied in Oschlies and Garc¸on (1999) and Oschlies et al. (2000). The traditional values for the linear and quadratic phytoplankton losses (F_p and F*_p) are taken from Oschlies (2001). No traditional parameter values exist for the temperature-dependent remineralization rates.

We, therefore, relate the previously published rates to the parameterization with tempera- ture by dividing the original rates by a factor of 3.5, which is the average value of the Eppley (1972) function over a temperature range from 5 through 28° C. The applied

traditional parameter set is listed in Table 6. The success of the optimization becomes evident when comparing the individual misŽ t contributions of the optimal Ž t and the traditional conŽ guration. For the total costs, the misŽ ts have been reduced by 25%. As desired, the simultaneous optimization has not only reduced the overall model data misŽ t but also the cost function contributions at each individual station. However, the minimum cost function still departs from the reference solution by a factor of 2.8. Such a failure might not seem very surprising, since the reference value could only be achieved under ideal conditions with a perfect model.

b. Optimal parameter values

A general feature of stochastic minimization algorithms is the very large number of cost function evaluations. Mapping all the cost function values generated during the minimiza- Table 5. Comparison of cost function contributions.The optimal conŽ guration resolves the contribu- tions of the best solution found after optimization. The reference values are derived from the optimal run with model results being extracted and processed similar to observational data at the dates of measurements. The traditional cost function relies on a model run with a formerly applied parameter set.

Cost function contributions

ConŽ guration J_{to ta l} J_{B A T S} J_{N A B E} JO W S - IN D I A

Optimal 60.59 17.71 25.04 17.84

Reference 21.65 4.72 14.17 2.75

Traditional 80.86 29.91 26.43 24.52

Table 6. The model’s parameter values. Traditional values refer to those typically applied in ecosystem models (see text). The best combination of parameter valuespˆ0is given together with approximated errorssˆ.

Parameters for optimization

n Symbolp_n Unit Traditional Optimal estimates:pˆ06 sˆ

1 m_m d²¹ 0.600 0.27060.033

2 a m²W²¹d²¹ 0.025 0.25660.036

3 F_m^p d²¹ 0.014 0.04060.013

4 F*_p m³mmol N²¹d²¹ 0.050 0.02560.014

5 k_N mmol Nm²³ 0.500 0.70060.010

6 k m²mmol N²¹ 0.030 0.04760.010

7 g d²¹ 2.000 1.57560.102

8 e m⁶mmol N²²d²¹ 1.000 1.60060.028

9 F_m^z d²¹ 0.009 0.01060.003

10 F*_z m³mmol N²¹d²¹ 0.200 0.34060.052

11 b 1 0.750 0.92560.060

12 g_m d²¹ 0.014 0.04860.020

13 w_s m d²¹ 5.00 18.0066.93

tion procedure can be used to gain some insight into the cost function’s shape. Figure 4 shows parameter guesses during the course of the optimization and their associated costs, projected onto the individual parameter axes. These Ž gures provide information on sensitivities of the cost function to parameter variations and hence on the relative importance of different parameters. Each subplot focuses on a single parameter that varies among the different forward integrations of the model and generates different cost function values (which will depend on values of the other parameters as well). Every combination of parameter values that occurred during the optimization is plotted together with the corresponding cost function value unless the costs exceed7_total595. The best parameter combination belongs to the lowest cost in the different subplots. Table 6 lists these best parameter estimates together with errors calculated according to Eq. 9. The faster the costs increase with distance from the best parameter value, the larger the cost function’s sensitivity with respect to the particular parameter. The subplots cover the full range of Figure 4. Parameter sets together with their lower limits in costs that occurred during the search process of all four optimizations (the three additional optimizations are included). Every abscissa shows the parameter’s value generated during the search process. The corresponding units can be gathered from Table 6. To the Ž rst subplot at the top all parameter projections are included (dots).

All other subplots simply show the lower limits of these projections. The fourteenth parameter is a control variable and has no effect on the model results.

parameter values between the upper and lower bounds. Under ideal conditions, in a sense that the optimization problem is well posed and the parameters are uncorrelated, all subplots would show sharp symmetric parabolas and the parameters would be fully constrained.

From Figure 4 it can be deduced how well each parameter is constrained. Among those parameters constrained rather well are the phytoplankton growth and loss parameters.

Estimates of the phytoplanktongrowth parameter (m_m) smaller than 0.27 d²¹cause a huge increase in costs. The same is true for estimates larger than the optimal value. The cost function sensitivity is less pronounced for the parameter of the photosynthetic efŽ ciency (a). For example, a value of 0.13 m²W²¹d²¹is signiŽ cantly lower than the best estimate of 0.256 m²W²¹d²¹, but it does hardly alter costs. Compared with commonly used values fora, the optimization value is larger by an order of magnitude. Possible reasons for this high estimate will be discussed in the next section. The well-constrained parameters describing the mortality of phytoplankton(FpandF*_p) also deviate from values which are conventionallyapplied: Higher optimal values are found forF_p, whereasF*_pis reduced by a factor of two. A relatively high half-saturation constant is obtained for the nutrient uptake (k_N5 0.7 mmol N m²³). The upper limit fork_Nis not reached but any value lower than 0.7 mmol N m²³ leads to an increase of the cost function. The relatively high optimal value of k_Nis somewhat surprising and we had initially expected much lower values because nutrient uptake in the NPZD-model re ects the uptake of both nitrate and ammonium, of which the latter seems to be relevant at low concentrations at the BATS and NABE sites during summer periods. We speculate that the optimal guess of kNis not independent of other parameters. The linkage between a few of the optimal phytoplankton parameter guesses, among which k_N is one, will be explained in more detail in the discussion section.

At Ž rst glance, the zooplankton parameters appear to be well constrained, but the two grazing parameters (g ande) and the assimilation efŽ ciency (b) are actually close to their respective upper bounds imposed to the optimization algorithm. The resulting estimates are close to those applied in other ecosystem models. To effectively constrain zooplankton parameters is difŽ cult with little zooplankton observations and no additional constraints on the grazing rates. In that case the prescribed limits of the parameter range become important. For example, the upper limit avoids the tendency toward an excessive grazing solution which would result in a modeled seasonal maximum of zooplankton biomass that is not resolved by the few available observations (e.g. Fennelet al.,2001; Schartauet al., 2001). In addition, any parameter combination which results in a strong reduction of zooplanktonbiomass enforces an increase in model-data misŽ t. Such a particular effect can be identiŽ ed for F_m^z. But there must also be a maximal tolerable zooplankton biomass, because the quadratic loss parameter for zooplankton (F*_z) reveals an increase of the costs as soon as its value is chosen to be lower than the best estimate. Producing a model solution with unreasonably high zooplankton biomass is likely to draw down the entire phytoplank-

ton biomass which leads to a data-model mismatch in chlorophyll and primary production after all.

The parameter for the remineralization of detritus (g) is hardly constrained by our cost function. However, the remineralization rate of detritus always depends on the best estimate of the sinking velocity (w_s) and vice versa. The optimal sinking velocity turned out to be signiŽ cantly higher than proposed for the traditional parameter conŽ guration.

Eventually, the remineralization rate increased by a similar factor as the sinking rate. This indicates, that the resulting remineralization proŽ le must be close to that generated by the traditional remineralization parameters.

The Ž nal “parameter” is a control which was introduced into the optimization process as well, but which has no effect on the cost function. With this control parameter we are able to determine whether the micro-genetic algorithm has converged to a particular solution for reasons other than being imposed by the shape of the cost function. The last subplot of Figure 4 illustrates that during the search process no preference toward a particular solution occurred. This is indicated by the range of Ž nal parameter guesses which still cover the entire search space within the prescribed upper and lower bounds, although a best combination of the other 13 parameters has already been identiŽ ed.

4. Discussion

a. Errors and sensitivity

Optimal parameter estimates are given by the minimum of the cost function. Neverthe- less, the form of the cost function and hence the location of its minimum depends on the actually available observations and their corresponding weights in the cost function. For a better interpretation of the optimal parameter estimates, information on precision, robust- ness and reliability of the optimization process is desired. To this end, a much simpliŽ ed Monte Carlo method (with Gaussian noise added to the observations and only 4 realiza- tions) has been applied. Following Eq. (9) this allows us to identify some uncertainties of the optimal parameter estimates which mainly re ect the robustness of the optimization algorithm. With this approach we can test whether we would retrieve a similar solution if our experiment was repeated with slightly different observational values (but of the same data type observed at the same dates).

Note that the additional (synthetic) solutions may be located quite close to the original point in the parameter space but the estimated errors may not recover the full error information which can be deduced from the cost function’s sensitivity to parameter variations. For example, the error estimate for the sinking velocity of detritus is close to 7 m d²¹for the optimal value of 18 m d²¹(Table 6), but Figure 4 reveals that the cost function hardly increases its value even for sinking velocity larger than 25 m d²¹. On the other hand, the estimated small error of the phytoplankton loss rate parameter Fmp is in accordance with the tremendous sensitivity of the cost function (see Fig. 4). As a consequence, the errors given in Table 6 must be jointly interpreted together with cost function sensitivities, as seen in Figure 4.

The above approach differs from another common error analysis based on the inverse of the Hessian matrix (Fennelet al.,2001; Vallino, 2000). The Hessian is the matrix of the second derivatives of the cost function with respect to the parameters, yielding the curvature or sensitivity. If computed at the minimum of the cost function, it can be used to estimate error bars for the optimal parameter estimates. A great advantage of computing the full Hessian matrix is that off-diagonal elements reveal correlations among the different parameters. In our case the dimension of the parameter space (513) would be small enough to allow direct computation of the full Hessian at little computational cost.

Unfortunately, when the model solution depends nonlinearly on the parameters in the vicinity of the cost function’s minimum, the inverse of the Hessian may not be a good approximation to the error covariance (Gunson and Malanotte-Rizzoli, 1996). In addition, the Hessian matrix can only be inverted as long as it remains well conditioned. To have a well-posed optimization problem is hardly the case when data are assimilated into an ecosystem model. In practice, though, this can always be satisŽ ed by adding ana priori constraint on the Ž rst parameter estimates to the cost function. However, we deliberately avoided sucha prioriassumptions on the parameter values.

b. Parameter estimates

The identiŽ cation of phytoplankton growth parameters is of great importance for algorithms which derive primary production rates from satellite data (Platt and Sathyen- dranath, 1999; Platt and Longhurst, 2000). The maximum growth rate of phytoplankton in theNPZD-model depends on the ambient temperature and on the optimal estimate of the growth rate parameter (m_m). Accounting for the typical temperatures at the three different sites, the resulting maximum growth rates within the uppermost layer are about 1.8 d²¹at the BATS location, 1.1 d²¹at the NABE site, and 0.8 d²¹at OWS-INDIA. The seasonal cycles of the modeled depth averaged maximum growth rates at the respective sites are shown in Figure 5. Evans (1999) used the same growth function but optimized a temperature independent maximum growth rate directly and obtained an estimate of 0.95 d²¹when assimilating only NABE-data into a FDM-model. This is in good agree- ment with our averaged rate for the NABE site. For the same location, the estimates of Fasham and Evans (2000) were higher, ranging from 2.02 d²¹ for the FDM-model to 1.27 d²¹for a modiŽ ed version with an additional diatom compartment. Nevertheless, the above comparison reveals that our growth rate estimate of 0.27 d²¹(times the temperature factor) is at the low end of estimates derived from data assimilation studies. Reasons for this low estimate of m_m must be discussed within a broader context, including other phytoplankton growth parameters as well, namely the initial slope parametera(photosyn- thetic efŽ ciency) and the half saturation constant for nitrogen uptake (k_N). For a similar ecosystem model, Fennelet al.(2001) demonstrated that the maximum growth rate can be negatively correlated with the initial-slope parameter (a) parameter and the half-saturation constant for nitrogen uptake (k_N). This is consistent with our results, where low growth rates coincide with high estimates of a and k_N. Our results of a are close to

0.25 m²W²¹d²¹, which is a factor of ten larger than the preferentially used value.

Estimates for a close to the traditional value were obtained by Fennelet al.(2001) for BATS data (a 5 0.0245 m²W²¹d²¹) and by Evans (1999) at the NABE site (0.035 m²W²¹d²¹). However, these two studies both comprised a priori parameter information in their cost functions which was a 5 0.025 m²W²¹d²¹. Other data- assimilation experiments found a tendency toward higher values ofaas well (e.g. 0.164 and 0.222 m²W²¹d²¹in Fasham and Evans (2000) or 0.93 m²W²¹d²¹in Hurtt and Armstrong (1999)). Similar values were estimated for aNPZ-model at the BATS location, ranging froma 50.173 to 0.688 m²W²¹d²¹(Schartauet al.,2001). The extremely high estimate of Hurtt and Armstrong (1999) was attributed to the circumstance that the initial time of the bloom at OWS-INDIA had to be Ž tted by their model correctly, getting a rapid bloom under low light conditions. Such a scenario is consistent with a time lag in simulated stratiŽ cation at the OWS-INDIA location (see accompanying paper Schartau and Oschlies, 2003).

To investigate which of the three locations enforced the high estimates ofa, Figure 6 resolves changes in the cost function that arise from varying two parameters,mmanda, respectively. Subplot 6A shows the shape of the cost function with the minimum at the optimal combination ofm_manda. The other subplots (6B–D) split up the cost function into the contributions that come from each location. The high estimates ofatogether with low growth rates result from the misŽ t contributions at the BATS and OWS-INDIA sites.

Furthermore, plot 6C reveals that a locally optimal parameter combination differs from the overall optimal solution (global minimum). This is of particular interest since it also proves that a single optimization at the NABE site would result in a different parameter set, e.g.

Figure 5. Modeled phytoplankton maximum growth rates averaged over the upper 35 meters.

Seasonal cycle at the BATS site (gray line), the location of NABE (black line), and at OWS-INDIA (dashed line).

estimates ofm_mandaclose to 0.4 d²¹and 0.12 m²W²¹d²¹, respectively. Ifais set to a traditional value, the cost function increases mainly due to the misŽ t contributions in DIN and Chla(Fig. 7). The increase in costs is dominated by changes in Chlaat the BATS site during all months, suggesting thata is constrained mainly by chlorophyll observations.

Naturally, it is the beginning of a phytoplanktonbloom which is sensitive to variations ina and the model is brought into better agreement with observations when the high estimate of ais applied. Low parameter values induce a signiŽ cant phase shift in the initialization of the phytoplankton bloom which generate large misŽ ts between model result and observa- tions (not shown). Apparently, the stratiŽ cation scenario at OWS-INDIA is not solely responsible for the optimal estimates, as discussed by Hurtt and Armstrong (1999). It is rather a general model deŽ ciency at the initial bloom phase. Here, a good example is found for a parameter estimate which will compensate for a model deŽ ciency. Apart from the initialization of the phytoplanktonblooms, high estimates ofaalso affect the estimation of the other growth parameters such ask_Nandmm. Sincearemains constant throughout the season, its high value causes light saturated growth conditions to great depths during the Figure 6. Two-dimensional variations of the initial slope parameter (a) and the growth rate parameter (m_m). Plot A shows the contours of the overall costs (as seen by the optimization algorithm). Subplots B–D reveal the individual contributionsto the cost function at every location.

Figure 7. Cost function contributions of DIN and Chlaat all three locations of interest. The misŽ t contributions are added up for each month at the respective site. Gray bars show the remaining misŽ ts after optimization. White bars indicate the misŽ ts which result when the initial slope parameter (a) is set to its traditionallyproposed value of 0.025 m²W²¹d²¹.

summer periods as well. But at depth between 50 to 100 meters, the modeled primary production rates are already higher than observed, as discussed in an accompanying paper Schartau and Oschlies (2003; this issue). The depth of the deep chlorophyll maximum during summer depends on estimates ofa, whereas the magnitude of the modeled deep chlorophyll concentration is largely controlled by the other growth rate parameters,k_Nand m_m, respectively. Hence, low growth rates are preferentially generated during the search process in order to minimize this individual misŽ t contribution.

High optimal estimates ofaresulted for models that compute a daily average of the light limited growth (Evans and Parslow, 1985) and do not account for diurnal variations in mixing. More precisely, such a daily average of the light-limited growth rate remains constant when integrating the phytoplankton equations throughout a single day (here 23 time steps for integration). Note that the net light availability for phytoplankton within the mixed layer can be enhanced when short-termed diurnal stratiŽ cation occurs while daylight becomes maximal. Such day and night changes in mixing depth can be signiŽ cant (Woods and Onken, 1982; Woods and Barkmann, 1986). When diurnal mixed-layer variations are neglected, the net light availabilityof cells trapped in the mixed layer at noon will be underestimated. We speculate that these effects are also present in our model, but are to some extent compensated by the parameter optimization. In order to receive improved estimates of a, we suggest a better accounting of the dynamics and physical- biological interactions during the initial bloom periods. Resolving diurnal mixing and a diurnal cycle for radiation is likely to make up for part of the models deŽ ciency.

Furthermore, if this model error is corrected, it would have a positive effect on other parameter estimates as well, such as onm_mandk_N.

The dynamics of the photosynthetic efŽ ciency itself may also be improved. For example, in Geider et al. (1998) a is modeled as a chlorophyll speciŽ c parameter, depending on the chlorophyll to carbon ratio. Furthermore, the authors suggest that the dynamical behavior of photosynthesis can be improved in the model when treatingaas the product of a light-absorption efŽ ciency and a maximum light-limited quantum efŽ ciency. The light-absorption efŽ ciency could then be treated as a function of the chlorophyll to nitrogen ratio while the quantum efŽ ciency becomes a function of DIN availability (it decreases under nitrate depleted conditions). Hence, low light and nutrient replete conditions would yield a high value forawhereas high light and low DIN concentrations would result in low values. Another sophisticated parameteriza- tion for the photosynthetic efŽ ciency has been suggested by Bissettet al.(1999), who have explicitly linked the pigment content with light absorption and actual energy utilization. Apart from the diurnal mixing effect, such improved dynamics is likely to account for a major fraction of our present model error.

While difŽ culties in Ž nding reliable estimates of the phytoplankton growth parameters can be attributed to model deŽ ciencies, the problems in estimating remineralization and export parameters are rather due to the constraints entering the cost function. In particular, the observational data are insufŽ cient to constrain the modeled phytoplankton sinks.

Phytoplankton losses can occur either via export of particle aggregates or by zooplankton grazing. Our estimates of the grazing parameters are close to the imposed upper bounds.

With these constraints we prevented the model to invoke a solution where high phytoplank- ton production rates and low chlorophyll concentrations are possible because of excessive grazing. However, the zooplankton parameters yield grazing rates between 1.5 and 1.6 d²¹ which may be too high because other processes responsible for a phytoplankton draw- down, such as cell coagulation and sinking of aggregates, are not sufŽ ciently accounted for in the model. If more zooplankton data were available, the choice of the upper bounds may become irrelevant and the modeled zooplankton biomass could be better optimized toward observations. With respect to the export, the two parameters of the detritus compartment are of special interest: the remineralization rate and the sinking velocity. The problem here is that both parameters cannot be estimated independently. Similar amounts of nutrients can be remineralized when sinking velocity and remineralization rate increase in the same proportions. Therefore, the estimated sinking velocity of 18 m d²¹needs to be considered in combination with the remineralization rate of 0.048 d²¹. It is noteworthy that the optimal estimates of the remineralization rate found here are constrained better than in an ecosystem model of the upper mixed layer only (Schartauet al.,2001). Doneyet al.(1996) proposed sinking rates of 10 m d²¹ together with a remineralization rate of 0.1 d²¹, matching a depth scale of remineralization of 100 m estimated from shallow sediment trap data by Lohrenzet al.(1992) for the BATS site. For the optimal parameters obtained here, the depth scale of remineralization becomes 375 m. The recycling of organic material of zooplankton is expressed as a linear temperature-dependent loss term in our model. When referring to the average temperature range at the three distinct locations, the  ux rates of nitrogen, from the zooplankton compartment back to its dissolved inorganic form, remain close to the traditionally temperature-independent constant rate. The best estimate for the linear remineralization rate of organic, presumably labile, nitrogen from phytoplankton to DIN is three times larger those values used in previous model studies at the respective sites.

Similar results were obtained by the data-assimilation study of Spitz et al.(2001). The authors pointed out the important role of bacteria for effectively utilizing dissolved organic nitrogen. Although this process is crudely resolved in theNPZD-model, our optimization result point toward high remineralization rates of labile organic material.

5. Summary and conclusion

A relatively simpleNPZD ecosystem model has been Ž tted to observations at three locations in the North Atlantic which re ect very different ecosystems. The minimum model-data misŽ t as measured by the minimum cost function value is about three times higher than could be reached by a model that was, apart from observational errors, fully consistent with the data. The reliability and robustness of the optimization procedure based on a genetic algorithm (mGA) is approved by additional optimizations with resampled data with added noise and by starting from different initial parameter guesses. For our

discussion of parameter values, we have found that the optimal combination of parameters, e.g. mm and a, at the NABE site can differ from values being optimal for the other locations. This is not a great surprise, but it clariŽ ed how the parameter estimates need to be interpreted. For example, if a different scaling approach was chosen and different weights were assigned, the misŽ ts at the NABE location could have altered the overall optimal estimates toward lower values for a and higher estimates for the growth rate.

Hence, the optimal parameter results are a compromise among all three locations, and the estimates are sensitive to the chosen scaling approach (which gives weights to the different locations). Nevertheless, one has to accept that ecosystem models always remain simpli-

Ž ed representations of the real biogenic environment. It will, therefore, hardly be possible to give reliable approximations of biogeochemical  uxes without carefully investigating the model’s parameter space as well. Such an exercise re ects the limits and weaknesses of the model assumptions. This study shows that data-assimilation experiments in ecosystem modeling are useful and can promote new model approaches. From our parameter optimization study we came up with the following conclusions:

a. The optimization resulted in unexpectedly high estimates of a. In fact, the optimized value of a is at the upper bound of the searchable parameter range. The model requires high values ofain order to reproduce observed chlorophyll concentra- tions mainly at the BATS site and, although less accentuated, at the location of OWS-INDIA. It is concluded that a high estimate a compensates for a model deŽ ciency. We suspect that the assumptions entering time-integrated functions for light limited growth might be inadequate.

b. High optimal values for the temperature-dependent phytoplankton loss rates are obtained. These high estimates are crucial for enhancing primary production rates by a rapid transformation of organic material back to its inorganic forms which are subse- quently available for phytoplankton growth.

c. The optimal parameter estimates must be interpreted as a compromise, albeit a reasonable one, among the three locations. This compromise will be sensitive to the chosen scaling approach, giving different weights to the individual locations. Hence, important parameters, such as those entering phytoplanktongrowth rate parameterizations, should be more dynamical in order to better account for the different ecosystem conditions found at the three different locations.

Acknowledgments.First of all, we thank everyone who supplied us with proŽ le measurements.

Many very helpful and constructive suggestions were made by the referees and we greatly acknowledge their efforts. We thank the staff of BBSR for making their very valuable data available through the US-JGOFS web-pages. We were generously supported by Roy Lowry and Polly Hadziabdic at the BODC. We thank Wolfgang Koeve, Joachim Hermann, and Mike Fasham who supported us with complementary data for 47N 20W and 59N 19W. The work presented in this paper was partly supported by the Deutsche Forschungsgemeinschaft (DFG).