Optimizing Models of the North Atlantic Spring Bloom Using Physical, Chemical and Bio-Optical Observations from a Lagrangian Float

The North Atlantic spring bloom is one of the main events that lead to carbon export to the deep ocean and drive oceanic uptake of CO 2 from the atmosphere. Here we use a suite of physical, bio-optical and chemical measurements made during the 2008 spring bloom to optimize and compare three different models of biological carbon export. The observations are from a Lagrangian float that operated south of Iceland from early April to late June, and were calibrated with ship-based measurements. The simplest model is representative of typical NPZD models used for the North Atlantic, while the most complex model explicitly includes diatoms and the formation of fast sinking diatom aggregates and cysts under silicate limitation. We carried out a variational optimization and error analysis for the biological parameters of all three models, and compared their ability to replicate the observations. The observations were sufficient to constrain most phytoplankton-related model parameters to accuracies of better than 15 %. However, the lack of zooplankton observations leads to large uncertainties in model parameters for grazing. The simulated vertical carbon flux at 100 m depth is similar between models and agrees well with available observations, but at 600 m the simulated flux is larger by a factor of 2.5 to 4.5 for the model with diatom aggregation. While none of the models can be formally rejected based on their misfit with the available observations, the model that includes export by diatom aggregation has a statistically significant better fit to the observations and more accurately represents the mechanisms and timing of carbon export based on observations not included in the optimization. Thus models that accurately simulate the upper 100 m do not necessarily accurately simulate export to deeper depths. Correspondence to: K. Fennel (katja.fennel@dal.ca)


Introduction
It is estimated that about 25 % of the global oceanic CO 2 uptake from the atmosphere takes place in the North Atlantic (Takahashi et al., 2009).A major contributor to this uptake is the North Atlantic spring bloom.According to classical theory (Sverdrup, 1953), the spring bloom is initiated when positive heat fluxes in spring cause a shallowing of the mixed layer allowing phytoplankton to remain long enough above the critical depth, where net primary production is positive, to induce positive growth.During the bloom phytoplankton typically grow rapidly, taking up nutrients near the surface.Export of Particulate Organic Carbon (POC) to the deep ocean results from sinking of organic particles primarily diatom aggregates and zooplankton fecal pellets.
The flux of POC into the deep ocean is dependent on the concentration and sinking rates of organic particles; the sinking rate, in turn, depends on the size and buoyancy of the particles.In the North Atlantic large diatoms often dominate the phytoplankton community at the beginning of the bloom (Sieracki et al., 1993) and constitute a major fraction of the sinking organic matter.Diatoms depend on dissolved silicate to build their frustules and take it up in approximately similar molar quantities as nitrate.Since silicate concentrations in the North Atlantic are lower than nitrate concentrations, silicate is typically the first nutrient to become depleted and to start limiting diatom growth (Allen et al., 2005).Physiological stress resulting from silicate limitation is known to increase the sinking rates of diatoms (Bienfang et al., 1982).By releasing extracellular polymeric carbohydrates, nutrientstressed diatoms increase the stickiness of their cell surfaces and thus support the formation of aggregates during collisions with other particles (Smetacek, 1985).The increase in size of these aggregates increases their sinking rates, which may exceed 100 m d −1 (Billett et al., 1983;Smetacek, 1985;Hegseth et al., 1995).As organic particles sink into the deep W. Bagniewski et al.: Optimizing models of the North Atlantic spring bloom ocean they are consumed by microbes releasing carbon and nutrients back into the seawater.The speed of sinking therefore determines at which depth the particles will be recycled.Diatoms are also thought to form dense cysts -a resting stage that is part of their life cycle -when silicate is exhausted (Smetacek, 1985).The cysts are capable of surviving for long periods of time in cold and dark deep waters (Malone, 1980).Depletion of silicate near the surface may prevent the growth of diatoms from late spring into late summer.During this time the phytoplankton community is dominated by smaller species and the microbial loop, which efficiently recycle nutrients and thus do not remove carbon efficiently from surface waters (Sieracki et al., 1993).
The North Atlantic spring bloom has been studied intensively, most notably during the JGOFS North Atlantic Bloom Experiment in [1989][1990] (Ducklow and Harris, 1993).Data collected during extensive cruises and from moored sediment traps significantly improved our understanding of the vertical carbon and nitrogen fluxes (Buesseler et al., 1992;Martin et al., 1993) and the plankton community's response to nutrient depletion (Sieracki et al., 1993).Modelling studies (Oschlies et al., 2000;Waniek, 2003) showed the need for more extensive and detailed observations for validating the models.Previous work also emphasized the need for a more complete spatial and temporal coverage of the bloom to grasp the variability caused by the observed mesoscale eddy stirring (Washburn et al., 1998;Mémery et al., 2005).These studies reveal the presence of small-scale phytoplankton patches, which create large differences in biogeochemical characteristics of surface water on horizontal scales as short as a few kilometers.Observations taken at fixed stations with moorings or during cruises do not follow the evolution of the bloom within a single water mass, but give merely a glimpse into the present state of every phytoplankton patch that passes through the site while transported by ocean currents.Satellite observations give a good snapshot of surface patchiness, but are scarce due to high cloudiness in this region and do not provide any information on subsurface distributions.
The North Atlantic Bloom 2008 (NAB08) experiment aimed at observing the bloom on the patch-scale using autonomous platforms.The experiment was conducted south of Iceland, near the 60 • N, 20 • W JGOFS site, using four seagliders and one Lagrangian float, each carrying a suite of physical, chemical and bio-optical sensors.The platforms were deployed in early April, well before the beginning of the bloom, and operated in the water until late June of 2008.The float remained in the mixed layer except for taking daily water column profiles to 250 m depth.This provided highresolution vertical data for one particular patch from prebloom conditions through full bloom development and decline.Ancillary measurements taken during four supporting cruises on the R/V Knorr and the R/V Bjarni Saemundsson provided physical, biological and chemical calibration data, information on the plankton community structure, and POC export measurements from sediment traps.Here the calibrated NAB08 float data are used for constraining and validating 1-dimensional ecosystem models.
Ecosystem models are mathematical representations of ecosystems in which biogeochemical processes are described with parameterizations (e.g., for the phytoplankton growth rate or mortality rate of zooplankton), based on our understanding of the system.The skill of a model can be defined by its ability to reproduce observations.Marine ecosystem models may vary greatly in complexity ranging from three to up to thirty and more biological state variables.Deciding how complex a model should be is perhaps the most difficult decision in model building (Hilborn and Mangel, 1997).According to Quine's theory of underdetermination, for every empirical data set there can be an infinite number of incompatible theories that explain it (Quine, 1975).While the more complex models may appear more realistic than the simple models, the increased complexity adds more parameterizations containing poorly known parameters (Denman, 2003).More complex models have more degrees of freedom, allowing to improve the fit with observations even when the data used in validation is sparse.However, as the number of parameters increases, it becomes more problematic to constrain these complex models with the available observations.Error analysis can help in determining which model parameters can be resolved given the available data.Studies by Matear (1995), Fennel et al. (2001) and Friedrichs et al. (2006) showed that many parameters in marine ecosystem models are highly correlated yielding large uncertainties of their optimized values.
Several marine ecosystem models have been developed in recent years for understanding, quantifying and predicting key biogeochemical processes in the oceans (e.g., Evans and Parslow, 1985;Fasham et al., 1990;Doney et al., 1996;Oschlies and Garc ¸on, 1999;Pondaven et al., 2000;Fennel et al., 2003a,b;Lima and Doney, 2004).Before routine use of parameter optimization techniques, models were typically optimized by tuning parameters manually until the model output "fits" the observations (e.g., Fasham et al., 1990).A preferable method is variational data assimilation, where model parameters are systematically perturbed to minimize a statistically based and objectively quantified misfit between the model and observations, thus providing an optimal parameter set (e.g., Matear, 1995;Spitz et al., 1998;Fennel et al., 2001;Friedrichs et al., 2006).
Here we set up a 1-dimensional physical model for the NAB08 location, based on the General Ocean Turbulence Model (GOTM), forced with atmospheric data, and nudged to temperature and salinity observations from the float in order to mimic the physical conditions along the float track.The model is coupled with three different biological models, which were optimized by assimilating the float data.The models were chosen to be similar in structure to models currently used for the North Atlantic.Our aim was to compare the three models in their ability to simulate the observed bloom and quantify the associated carbon export, and to investigate the importance of including diatom aggregation triggered by low silicate concentrations.

In situ observations
Data used for validation and assimilation were collected by a heavily-instrumented Lagrangian mixed-layer float (D'Asaro, 2002), which operated at the NAB08 site (Fig. 1) from Year Day (YD) 95 (4 April) to YD 146 (25 May).The float took one vertical profile daily from the surface to 250 m; otherwise it drifted passively in the mixed layer.The properties measured by the float are given in Table 1.Ancillary measurements for in situ calibration of the float sensors were made during a cruise on the R/V Knorr from YD 123 to 142 (2-21 May) and the platform deployment and recovery cruises on the R/V Bjarni Saemundsson from YD 94 to 96 and YD 178 to 181 (3-5 April and 26-29 June).The measurements include extracted chlorophyll, POC, oxygen, nitrate and silicate from bottles, as well as CTD measurements of oxygen, chlorophyll fluorescence and beam attenuation.Chlorophyll fluorescence and beam attenuation from the float were used as proxies of chlorophyll and POC concentrations, calibrated by the it in situ measurements.Details of the measurements, associated errors and the calibrated data can be found in calibration reports and archives at http://osprey.bco-dmo.org.The calibrated float data were averaged in bins of 1 m vertically in the top 100 m and 1 day temporally.Silicate measurements collected within 10 km of the float were used.Bins without data points were filled by interpolation (except for silicate).
The estimated errors in these data products are listed in Table 1.Two error sources are considered.The sampling er-ror measures the uncertainty due to measuring N points in a bin rather than all possible points.It is computed from average standard deviation in a bin, divided by the square root of the average of 1/N, a value always between 0.8 and 1. Measurement error is due to calibration uncertainty of the sensor and interpretation uncertainty of the proxy measurements.This error includes both random and bias components.The first is included in the sampling error, but the second is not.To be conservative, these two error sources are considered independent and their squares summed to get the total error.The measurement error dominates for chlorophyll, nitrate and POC and approximately equals the sampling error for oxygen.For silicate, the error is estimated from the average standard deviation of the bottle data values at the sampled depths relative to a multi-day smoothed spline fit to these data.
In addition to the data used directly in model initialization and assimilation (listed in Table 1), other measurements from the experiment gave additional insight and were used to inform our choices of model structure and as it a priori knowledge in the data assimilation procedure.Specifically, spikes in optical measurements appeared as silicate became depleted and moved vertically over the course of a few days (Briggs, 2010) indicating the formation of sinking diatom aggregates.Samples captured by PELAGRA floating sediment traps contained large numbers of Chaetoceros (diatom) chains and resting cysts (Martin et al., 2011) and were estimated to sink at about 50 to 100 m d −1 (Briggs, 2010).This export event was accompanied by a rapid shift in the plankton community composition, which was dominated by diatoms before the export event and dominated by picoeukaryotic phytoplankton after the export event (M.Sieracki, personal communication, 2010).

Model description
The biological models used in this study are embedded into a 1-dimensional physical model (the General Ocean Turbulence Model or GOTM; Burchard et al., 1999).The physical model describes the top 200 m of the ocean with a vertical resolution of 1 m, and uses the k-epsilon mixing scheme.The physical model is forced with wind speeds, air pressure, air temperature and humidity from the NCEP/NCAR reanalysis data set (Kalnay et al., 1996) except for YD 123 to 142 (2-21 May) when wind speeds from the R/V Knorr's meteorological tower were used.Solar radiation was calculated from the float's PAR sensor by extrapolating to the surface using a daily varying attenuation coefficient.Since horizontal transport processes are not resolved in the model, the model temperatures and salinities were nudged to their corresponding observations.A nudging time scale of 6 h resulted in close agreement between model-simulated and observed mixed layer depths.As a result, the model mimics the physical conditions along the track of the Lagrangian float.Silicate was the only nutrient that reached limiting levels during our simulation period.Nitrate does not reach limiting concentrations, but was measured and is therefore included in the model.Iron was not included because it is unlikely to be deficient in the study area in the spring (Martin et al., 1993;Fung et al., 2000), although iron may limit phytoplankton growth in the summer (Nielsdóttir et al., 2009).
We compare three biological model variants.The first (referred to as 1p1s) includes one phytoplankton group, Phy, two nutrient groups (dissolved inorganic nitrogen, DIN, and silicate, Si), two types of detritus (detrital nitrogen, Det N , and detrital silica, Det Si ) and zooplankton, Zoo.The other two variants (referred to as 2p1s and 2p2s) include two phytoplankton groups, representing diatoms, Dia, and small phytoplankton, Phy.The 2p2s model also includes diatom aggregates and cysts, Cys, with a different sinking speed than diatoms (Fig. 2).Biological model variables and parameters are listed in Table 2.The simplest model (1p1s) is similar in structure to the model of Oschlies and Garc ¸on (1999), which has been used extensively for the North Atlantic in many follow-up studies (e.g., Oschlies et al., 2000;Oschlies, 2002), although here we included chlorophyll as a separate state variable to allow for photoacclimation and included the uptake of silicate by phytoplankton in addition to the uptake of nitrate.The other two model variants (2p1s and 2p2s) are similar in structure to the model of Lima and Doney (2004), which has also been applied to the North Atlantic, in that diatoms are represented as a separate phytoplankton group.
Particulate organic nitrogen is calculated as the sum of small phytoplankton, diatoms, cysts, zooplankton and detrital nitrogen and is compared with observed POC assuming the Redfield ratio.Sinking organic matter leaves the model domain upon reaching the bottom boundary.The sources and sinks of the biological variables are given below.The model was run from YD 111 (20 April), before the observed start of the bloom, to YD 145 (24 May), when the float was recovered.
where the rates µ, m, g and w represent growth, mortality, grazing and sinking rates, respectively.Phytoplankton mortality is assumed to represent all non-predatory death, e.g. by viral cell lysis or physiological senescence (Franklin et al., 2006).The phytoplankton growth rate depends on nutrient concentrations, water temperature, T , and the amount of photosynthetically active radiation, PAR, according to   oxygen-to-DIN stoichiometry of organic matter production and respiration zooplankton mortality rate R dead -initial dead organic matter per total organic matter R phy -initial phytoplankton per living organic matter R dia -initial diatoms per all phytoplankton Here µ max is the temperature-dependent maximum growth rate µ max (T ) = 1.066T µ 0 with µ 0 as the maximum growth rate at T = 0 • C (Eppley, 1972).Nutrient limitation is represented by the Michaelis-Menten parameterization with k N and k Si as the half-saturation concentrations.Lightlimitation is parameterized as given in Evans and Parslow (1985) and has the form where P refers to Phy and Dia, α is the initial slope of the photosynthesis-irradiance curve and PAR is the photosynthetically active radiation.PAR decreases exponentially with depth, z, according to Here I 0 represents the total incoming radiation just below ocean surface and φ = 0.43 is the fraction of I 0 that www.biogeosciences.net/8/1291/2011/Biogeosciences, 8, 1291-1307, 2011 lies in the photosynthetically active spectral range (Kirk, 1994).The coefficients K w = 0.059 m −1 and K Chl = 0.041 (mg chl) −1 m 2 represent light attenuation due to water and chlorophyll, respectively.The biological dynamics of zooplankton is given by where g is the grazing rate, β is the assimilation efficiency (the term (1 − β)g accounts for waste produced from undigested food, which forms fecal pellets and enters detritus), ε 1 is the rate of ammonium production through excretion, and ε 2 is the mortality rate.Zooplankton grazing is parameterized using a modified Ivlev equation (Franks et al., 1986) and depends on small phytoplankton concentration with g max as the maximum grazing rate and λ as the Ivlev grazing coefficient.
The nutrients (DIN and Si) are produced through detritus remineralization, r, and zooplankton excretion, and are removed through phytoplankton uptake according to and Here R Si:N is the Si-to-N stoichiometry of diatoms.Detrital nitrogen is formed from dead phytoplankton and zooplankton, and from zooplankton fecal pellets.Detrital silicate is formed from dead diatoms only.Detritus is remineralized at the rates r N and r Si for detrital nitrogen and detrital silicate, respectively.Both detrital groups are fractions of the same detrital pool and, thus, sink at the same rate w Det .
Chlorophyll of small phytoplankton and diatoms (Chl 1 and Chl 2 , respectively) is affected by the same processes as phytoplankton, except that its production is controlled by the photoacclimation factor ρ Chl .
The photoacclimation factor accounts for the variation in the chlorophyll-to-phytoplankton biomass ratio with light availability (Fennel and Boss, 2003).It causes increased chlorophyll production under low light conditions and is determined following Geider et al. (1997) as where P represents Phy and Dia, and θ max is the maximum chlorophyll-to-nitrogen ratio.
Oxygen is produced during photosynthesis and consumed by zooplankton metabolism and detritus remineralization.Therefore, oxygen production is directly proportional to changes in DIN, scaled by the oxygen-to-DIN stoichiometry (R Oxy:DIN ) 16) At the ocean surface, oxygen concentrations are affected by gas exchange with the atmosphere, F airsea , which is parameterized following Woolf and Thorpe (1991) as Here K is the coefficient for oxygen exchange, O sat is the temperature and salinity dependent oxygen saturation concentration calculated based on Garcia and Gordon (1992), U 10 is the wind speed at 10 m above the sea surface and U 10x is the wind speed at which the equilibrium oxygen supersaturation is 1 %.Since F airsea has a 30 % uncertainty, we included the scaling parameter k Oxy , initially set to 1, which is included in the optimization described below.After Nightingale et al. (2000) the oxygen exchange coefficient K is where S C is the Schmidt number calculated following Wanninkhof (1992).
The variable Cys represents rapidly sinking diatom aggregates that include diatom chains, broken frustules and cysts.They are formed from diatoms at time t Cys (the time when silicate concentrations in the surface layer fall below a critical level, Si Cys ).All diatoms are instantly transferred into Cys, which sinks at the rate w Cys and does not decay due to grazing, mortality or remineralization.Cysts are assumed to preserve the chlorophyll from diatoms, hence Chl 2 also sinks at that rate.This treatment of diatom sinking upon depletion of silicate is similar to that of Pondaven et al. (2000).The Dirac delta function, δ, is used here to indicated the transfer of diatoms into cysts at time t Cys .
The model was initialized with observed values of chlorophyll, Particulate Organic Nitrogen (PON), nitrate, silicate and oxygen.Since no observations were available for phytoplankton, zooplankton or detritus concentrations, initial PON values were distributed between these variables using the adjustable parameters R dead , R phy and R dia as follows The values for R dead , R phy and R dia were optimized as described below.

Variational data assimilation
We use variational data assimilation to optimize the biological model parameters and compare the performance of the optimized model variants in reproducing our observations.As in non-linear least squares fitting, the squared misfit between the observations and their model counterparts is measured by a cost function which is to be minimized.The cost function has the form where M is the number of data types (here We used an additional term, F R , in the cost function of the 2p1s and 2p2s models in order to consider the qualitative knowledge of phytoplankton species composition.Before YD 133 (12 May) the term F R contributes to the cost function if diatom biomass is less than twice the biomass of small phytoplankton.After YD 137 (16 May) F R contributes to the cost if diatom biomass is more than half of small phytoplankton biomass.Specifically, where R spec is defined as The cost function used for the optimization of the 2p1s and 2p2s models becomes F + W R F R where W R is a dimensionless weight and was chosen such that the cost contributed from the term W R F R is roughly equal to the contributions of each of the other variable types.
The cost function was minimized using the gradient descent routine by Gilbert and Lemaréchal (1989).This routine requires the gradient of the cost function as input, here approximated numerically by perturbing every parameter p by a small p of 10 −5 .Since gradient descent methods can not distinguish local from global minima, we applied the minimization repeatedly from randomly chosen initial parameters to ensure that the global minimum of the cost function is found.The initial parameters were drawn from uniform distributions bounded by the minimum and maximum parameter values given in Table 3.For the 1p1s and 2p1s models, the minimization routine was run 50 times, with all parameters varied simultaneously.About half of the minimizations appear to have stopped at local minima while the rest converged on one absolute minimum which we consider the global minimum.The corresponding parameter values agree to within ∼10 −4 of their absolute value.Due to higher complexity and possibly a much greater number of local minima caused by the sharp "switch" from diatoms to cysts, the 2p2s minimization problem is poorly conditioned and not all parameters could be varied at once.For this model subsets of ∼10 parameters resulted in a well conditioned minimization problem and were optimized simultaneously.By cycling through different parameter subsets we optimized all parameter values and gradually reduced the cost function.To reach a convincing global minimum, the 2p2s model was minimized ∼1500 times for each parameter set.
Two biological parameters, the sinking rate of diatom cysts, w Cys , and the half-saturation concentration for DIN uptake, k N , were not optimized.w Cys was set to −55 m d −1 , based on backscatter observations from the seagliders as described in Briggs (2010).k N was poorly constrained by the observations since observed DIN concentrations were always well above typical values for k N .
We estimated the uncertainty in the optimal cost function values due to observational errors using a Monte Carlo analysis.Specifically, we perturbed the observations by adding random noise to each data point drawn from normal distributions with a mean of zero and a standard deviation equal to the standard error for the corresponding variable type (see Table 1).We created 1000 such randomly perturbed data sets and calculated the resulting cost function values for each of the three optimal models.The standard deviations of these cost function values provide estimates of the standard error of the optimal costs for the three models and are given in Table 4 for the total cost as well as for the cost contributions from each variable type.

Hessian matrix
In biological models parameters are typically correlated with each other, making certain combinations of these parameters poorly determinable by the data.Also, there may be little or no information about some of the parameters.The Hessian matrix of the cost function can be analyzed to estimate these correlations as well as a posteriori parameter uncertainties and the conditioning of the minimization problem (Thacker, 1989;Fennel et al., 2001).We approximated the Hessian of the cost function numerically by perturbing the model parameters and calculating the second derivative of the cost function for each perturbation.For the 2p1s and 2p2s the cost function including the F R -term was used.
Near the global minimum, the inverse of the Hessian matrix provides a good approximation of the error covariance matrix for the independent model parameters (Thacker, 1989).The condition number of the Hessian matrix, calculated as the ratio of the matrix's largest to smallest eigenvalue, shows how singular the minimization problem is.In the case of a large condition number, corresponding to an ill-conditioned matrix, the convergence rate of the minimization algorithm is slow.For the 2p2s model we optimized only subsets of parameters with a condition number smaller than 2 × 10 4 .The a posteriori errors of model parameters are the diagonal elements of the inverse of the Hessian matrix.

Optimization
The variational assimilation method described in the previous section was applied to the three model variants and resulted in the following minima of the cost function F : 5.6 ± 0.12 for the 1p1s model, 5.0 ± 0.12 for the 2p1s model, and 4.6 ± 0.11 for the 2p2s model (the contribution from term F R is excluded here as it is not used for the 1p1s model, see Table 4).In other words, the best fit between model and observations was found for the 2p2s model.Also, the fit for the 2p1s model is better than for the 1p1s model.The value of F + F R is 5.6 and 5.3 for the 2p1s and 2p2s models, respectively, again indicating that the 2p2s model gives the best fit.
The optimized biological model parameters are shown in Table 3.Most of the parameters are very similar between the 2p1s and 2p2s models.There is generally a better agreement for parameters with lower a posteriori errors (e.g.µ 0 , α, β, r dia ) than for parameters with higher errors (e.g.w Phy , k Si , λ), confirming that their values are indeed more tightly constrained.Large mismatches are found for the parameters directly affecting zooplankton (λ, g max , ε 1 , ε 2 , r phy ).For the 1p1s model some parameters differ from those of the 2p models, including some directly related to phytoplankton (m Phy , α Phy ).The only parameter measured during the NAB08 experiment was α.Productivity measurements during the Knorr cruise found a mean value of α = 0.034 mmol N mg −1 Chl a d W m −2 (K.Gudmundsson, personal communication, 2010), which represents the community mean.

Model results
Model temperature (Fig. 3) and salinity follow the observations closely, because of strong nudging.As a result, the model-simulated and observed mixed layer depths (MLD) are highly correlated with an R 2 of 0.89, where the MLD is defined as the depth at which the density difference to the sea surface is 0.01 σ θ or higher.A brief, first shallowing of the MLD to less than 20 m took place during the second half of April, but the MLD deepened soon afterward to ∼100 m and remained deep until early May (YD 127).For about 10 days from early to mid-May (YD 127 to 138) the MLD was shallow at or above 30 m, then deepened to about 30-40 m thereafter.Cooling and freshening of the water below 50 m occurred after YD 130 (9 May) and were likely caused by horizontal advection of a different water mass.
Results from the optimized biological models are shown with observations in Figs.4-8.Observed chlorophyll concentrations are initially at ∼0.5 mg Chl m −3 (YD 111) and increase with the bloom until YD 133 (12 May).During these 23 days surface concentrations increase 7-fold and reach a maximum of ∼ 3.5 mg Chl a m −3 (Fig. 4).The growth is interrupted around YD 120 (29 April), during a storm which deepened the mixed layer (see Fig. 3) decreasing surface concentrations and mixing chlorophyll to greater depths.A rapid decline, lasting ∼3 days, follows and depresses chlorophyll further to values only slightly above pre-storm values.The model-predicted chlorophyll evolution has similar patterns to the observed chlorophyll, although surface concentrations are overestimated initially during the bloom, especially in the 1p1s and 2p1s models.The 2p2s model matches the observed chlorophyll best, and predicts a rapid decrease in chlorophyll after the bloom at all depths matching the observations closely.The 1p1s and 2p1s models predict a more gentle decrease.
PON follows a similar pattern as chlorophyll (Fig. 5).During the bloom observed PON concentrations grow from ∼ 0.7 mmol N m −3 to ∼ 3 mmol N m −3 at the surface, and decline rapidly after the bloom to ∼ 1.5 mmol N m −3 .The model-predicted PON evolution follows the observations The observed DIN initially has concentrations of 13 mmol N m −3 .As it is consumed by phytoplankton its surface concentrations decrease to ∼ 8 mmol N m −3 in 23 days during the bloom (Fig. 6).Surface DIN briefly increases around YD 120 and 130, due to deepening of the mixed layer.At the end of the bloom the surface concentrations stabilize and slowly increase upon mixing with deeper layers.The three models match the observed DIN very closely during the bloom, but start to diverge afterward.In the 1p1s model the bloom ends too early, while in the 2p1s and 2p2s models the DIN concentrations continue decreasing for a few days  7).Below 50 m the oxygen increase due to photosynthesis and mixing with upper layers is counterbalanced by respiration.Again the models match the observations closely until the end of the bloom when their predictions start to diverge.And again, the 2p2s model gives a better match than the other two models.
Silicate concentrations were not measured autonomously and are available only during the R/V Knorr cruise, i.e. from

Plankton composition
In both 2p models diatoms dominate biomass during the bloom (Fig. 9).The 2p2s model predicts larger concentrations of both phytoplankton groups than the 2p1s model, until all its diatoms sink out.In the 2p1s and 2p2s models, diatom concentrations start declining around YD 130 as silicate limitation begins to hamper diatom growth.Diatoms are completely removed from the 2p2s model on YD 133 when they are converted into aggregates and cysts and sink.In the 2p1s model the decline is more gradual.The models' small phytoplankton concentrations peak at the time of diatom decline and decrease subsequently after YD 134, mostly due to increased grazing pressure.
In the 2p1s and 2p2s models, zooplankton concentrations decrease during the first half of the simulation, when small phytoplankton concentrations are still very low, which seems unrealistic.In the 1p1s model, grazing begins earlier and increases more gradually than in the other models.Although no zooplankton data were assimilated and the models use very different zooplankton parameter values, zooplankton biomass is very similar in all three models.
Detailed quantitative information about the plankton community structure is not available, but some 10 m bottle samples were taken during the R/V Knorr cruise and analyzed by flow cytometry for picoplankton and by imaging-in-flow (FlowCAM, Sieracki et al., 1998) for microplankton.The ratio of small phytoplankton to total phytoplankton carbon increased from 0.2-0.4 on YDs 126-128 (5-7 May) to ∼1 after YD 137 (16 May) in these measurements (M.Sieracki, personal communication, 2010).In the 2p1s model the ratio increased from ∼0.2 on day 125 (4 May) to 0.7 by day 134 (13 May) and remained at 0.7 thereafter.In the 2p2s model the ratio was similar to the 2p1s model up to YD 134, but then increased to 1 and was in better agreement with the observed data than the 2p1s model.

POC fluxes
POC fluxes were calculated directly from the concentrations and sinking rates of model variables, using the Redfield C:N ratio.The estimated POC export at 100 m was similar for the three model variants ranging from ∼1.35 mol C m −2 in the 1p1s model to ∼ 1.61 mol C m −2 in the 2p2s model (integrated for the simulation period from YD 111 to 145; Table 5).At 100 m depth diatom aggregates and cysts represent 32 % of the export in the 2p2s model, while diatoms represent 21 % of the export in both 2p models (Table 5).Detritus contributes 94 %, 78 % and 47 % in the 1p1s, 2p1s and 2p2s models, respectively.The large contribution of detritus to sinking organic matter is mostly due to its high sinking rates.Export related to diatoms is most important in the 2p2s model due to the switch to fast sinking aggregates and cysts.POC export is relatively constant throughout the experiment, except for the major event of sinking diatom aggregates and cysts in the 2p2s model on YD 134 (Fig. 10).Before this event the 1p1s and 2p2s models predicted very similar POC export, despite having different compositions of the sinking matter (Fig. 10).No measurements of the sinking flux were used in data assimilation, hence, the simulated POC export is informed only by its imprint on PON and DIN, which were assimilated.
In order to estimate export fluxes at greater depth we extended the optimal models vertically from 200 m to 600 m depth.The POC fluxes at 600 m were much lower than at 100 m and differed markedly between the three models (Fig. 10).At 600 m export predicted by the 2p2s model was 2.5 times larger than in the 1p1s model and 4.5 times larger than in the 2p1s model.80.2 % of the export in the 2p2s model is from sinking aggregates and cysts, which reach 600 m depth during the last days of the simulation resulting

Discussion
We optimized three different variants of a 1-D ecosystem model for the North Atlantic spring bloom using a suite of high-resolution observations from a Lagrangian float.The biological models are very sensitive to changes in vertical mixing, as every deepening of the mixed layer dilutes and redistributes biomass and nutrients.We achieve a realistic representation of the mixed layer evolution by nudging strongly to the observed temperature and salinity fields.The model essentially replicates the observed fields of these variables.This alleviates one of the main limitations of the 1-D set up, namely the lack of horizontal buoyancy fluxes, such as the observed cooling and freshening from below in late May (see Fig. 3, near YD 139), although nudging does not account for horizontal advection and mixing of biomass and nutrients.One potential outcome of parameter optimization studies is for the optimization to fail because the model does not adequately represent the system that is studied, or because the available observations are not adequate to constrain the model or both.This has been the case in several previous biological optimization studies.For example, Fennel et al. (2001) and Schartau et al. (2001) optimized simple mixed layer models for the Bermuda Atlantic Time Series station and concluded that their models were not adequate to fit the observations.Schartau et al. (2001) reported ecologically unreasonable optimal parameter values, while Fennel et al. (2001) found large parameter uncertainties.Ward et al. (2010) optimized 10 parameters for two models of differing complexity for the Arabian Sea and found that the available observations did not contain enough independent information to constrain all 10 parameters.Ward et al. (2010) concluded that "[...] even simple marine biogeochemical models are currently underdetermined by observations at oceanic time-series sites [...]."The optimizations described here were successful in that all biological parameters were optimized and had plausible values (many with small uncertainties) and that all models reproduced the main features of the bloom, and only varied in some details (e.g.maximum chlorophyll concentration, rate of the bloom decline, or amount and composition of carbon exported).
Model performance can be quantified by the squared mismatch between observations and their model counterparts (i.e. the cost function, Eq. 24).By this measure the most complex model (2p2s), which includes the formation of fast sinking diatom aggregates under silicate limitation, shows the best performance (i.e. it has the smallest cost or misfit).The fit between optimized model and observations improved with increasing complexity (by about 10 % by adding a second phytoplankton functional group to the 1p1s model, and by another 8 % when diatoms were allowed to form fast sinking aggregates and cysts in the 2p2s model) as one would expect.While the better fit for the more complex model could simply result from the added degrees of freedom, the most straightforward interpretation, in our case, is that the more complex model fits the data better because it better represents the biological processes.Other data, not used in the optimization, supports this interpretation as is described below.It would have been desirable to test the optimized models against unassimilated data.However, this was not possible using NAB08 observations due to the lack of independent data.Two Lagrangian floats were deployed during NAB08, which would have resulted in two independent trajectories, but one float malfunctioned early in the experiment.It has been suggested that any parameter optimization study should include an error analysis that estimates the uncertainty of the optimal parameters (e.g., Matear, 1995;Fennel et al., 2001).However, a posteriori errors are frequently not reported in biological parameter optimization studies.A by no means complete list of exceptions includes the studies by Matear (1995), Prunet et al. (1996), Fennel et al. (2001) and Faugeras et al. (2003).Here, 12 out of 25 biological parameters had small uncertainties (<15 %) for the 2p1s model and 7 out of 19 for the 1p1s model.In contrast, Matear (1995) was able to determine between 3 and 5 biological parameters with uncertainties of <15 % when optimizing three different models for Station P in the North Pacific, and Fennel et al. (2001) found uncertainties larger than 15 % for all parameters when optimizing a model for the Bermuda Atlantic Time Series.Here the parameters directly affecting phytoplankton (e.g.mortality m, maximum phytoplankton growth rate µ 0 , initial slope of the P − I curve α) were typically more certain than the parameters directly affecting zooplankton (e.g.maximum grazing rate g max , Ivlev grazing coefficient λ, excretion ε 1 and r phy , which determines the initial concentrations).This is consistent with previous studies.A large uncertainty implies that the corresponding parameter is not constrained by the observations and that the cost function is not sensitive to changes in the parameter.Parameters related to phytoplankton growth are well determined (most of the optimized values were similar for the different models and had errors of less than 15 %), because phytoplankton was observed with two proxies (PON and chlorophyll) and had the largest dynamic range of all biological groups, leaving a distinct trace in nutrients and oxygen.Similarly, oxygen parameters (r Oxy:DIN , k Oxy ) were found with little uncertainty.Parameters related to zooplankton had larger uncertainties, even though the simulated zooplankton biomass was very similar between the models, with low concentrations during and a rapid increase near the end of the diatom bloom and stable concentrations thereafter.The agreement of zooplankton concentrations between models after YD 137 (16 May), when the system appears to have been in a steady state, suggests that zooplankton biomass was constrained by the observations at this time.However, the optimized values of the parameters related to zooplankton (e.g., λ, g max ) were very different between the models and had large uncertainties.Hence, even if zooplankton biomass can be constrained during periods of relatively steady state this does not imply that rates are constrained as well.The parameters k Si and m Dia (m Phy in 1p1s model) determine the timing and the steepness of the diatom decline and differ in all models.Without the fast sinking of diatoms, the 1p1s and 2p1s models require higher diatom mortality to simulate the observed decline of chlorophyll and PON at the end of the diatom bloom, and a lower half-saturation concentration of silicate to ensure that this decline does not happen too early.
On the other hand, w Phy is very small and had only a minor effect on the model, therefore it is associated with a large error.
Although all three model variants provided a satisfactory fit with the data in a least squares sense, the 2p2s model is superior in reproducing several qualitative aspects of the observations.Observations of community structure by Flow Cytometry and FlowCAM analyzes, show a rapid shift from bloom conditions, dominated by diatoms, to post-bloom conditions dominated by smaller cells (M.Sieracki, personal communication, 2010).Only the 2p2s model reproduces this rapid change.Observations also clearly show that the disappearance of the diatoms occurred at the time of silicate exhaustion through aggregation and rapid sinking similar to that described by Smetacek (1985), Bienfang et al. (1982), Sieracki et al. (1993) and Dale et al. (1999).The sinking aggregates were captured in sediment traps (Martin et al., 2011) and shown to be rich in diatom aggregates and cysts.The distribution of the aggregates was mapped using spikes in optical measurements of backscatter, beam attenuation and chlorophyll fluorescence (Briggs, 2010) and shown to be widespread at the occurrence of silicate depletion.Only the 2p2s model reproduces this observed rapid export of diatoms.
The three model variants predict similar POC export at 100 m depth, because all three are tuned to match the observed net community production above 100 m depth and structured so that this production must be exported.The value of 1.6 mol C m −2 over 34 days, or 47 mmol C m −2 d −1 , in the 2p2s model lies within the estimated range of 22.5-67 mmol C m −2 d −1 based on a Th-234 analysis of sediment trap data performed during the Knorr cruise (Martin et al., 2011).The simulated export values are in reasonable agreement with values obtained during the JGOFS North Atlantic Bloom Experiment by Martin et al. (1993) of 39 mmol C m −2 d −1 , by Boyd and Newton (1995)  Although the POC export at 100 m is necessarily similar, differences between the three model variants produce large differences in the mechanisms of export and thus in its timing and depth.In particular, export in the 1p1s and 2p1s models occurs primarily though slowly sinking detritus, so that significant remineralization occurs.In contrast, the rapid sinking of aggregates in the 2p2s model reduces the amount of remineralization and thus increases the efficiency of POC export to depth.Accordingly, POC export at 600 m is much larger for 2p2s than for the other two models.Note that the depth of wintertime mixing in this region always exceeds 100 m, so that carbon exported to 100 m is easily mixed back to the surface during the following winter.Mixing to 600 m is much less common, so that export to depth is more likely to result in sequestration of this carbon for many years.
Thus in this study, three models with small differences in formal least squares error have large differences in their qualitative reproduction of export mechanisms and in their predicted export rates at 600 m.This highlights the difficulties in finding, or even defining, a "best" model by formal methods alone.Although the use of high quality data in the upper 100 m results in highly accurate estimates of many model parameters, the optimization scheme used here does not penalize models that fail to reproduce the observed rapid export.Additions to the model cost function that included such penalties could easily be added in a manner similar to Eq. ( 25) and, with arbitrary weight, could induce an arbitrarily large impact on the accuracy of the model depending on the importance that the modeler places on reproducing this feature.A formal definition of "best" without such a penalty is equivalent to giving it no weight.
More importantly, this study demonstrates that models that accurately simulate the upper 100 m do not necessarily simulate export to deeper depths accurately.It appears that detailed data extending through the twilight zone will be necessary to develop models that accurately predict export through this zone.

Conclusions
Using a suite of high-resolution physical, bio-optical and chemical observations from the NAB08 Lagrangian float, we were able to constrain the evolution of vertical mixing in a 1dimensional ocean model and to optimize three different biological model variants, after supplementing the observations with ship-based measurements of silicate.Since inorganic nitrogen never reached limiting concentrations, the models' response to limitation by silicate, which can not be measured autonomously, was crucial for capturing the termination of the spring bloom.In contrast to previous biological optimization studies, many of the optimal parameters had small uncertainties (less then 15 %), which suggests that biological model development and calibration would greatly benefit from a more widespread availability of bio-optical and chemical float observations.Consistent with earlier studies, the most uncertain parameters were related to zooplankton and vertical sinking.When comparing the three biological models, the fit between models and observations improved slightly with increasing complexity, but the most complex model (2p2s) was also the most difficult to optimize (i.e.not all parameters were optimized simultaneously).Given the only slight differences in misfit between the models and observations for the three variants, none can be justifiably rejected based on misfit alone.However, only the model variant that includes diatom aggregate formation and sinking qualitatively simulates the observed collapse of the diatom bloom by this mechanism and the resulting observed rapid export of carbon to depth.Thus, increased complexity in biological models may still be justified based on ecologi-cal or theoretical considerations.For example, if the models described here were used in the context of a basin scale biogeochemical circulation model, the variant with rapid vertical sinking of diatom aggregates would likely lead to different conclusions about the efficiency of biological carbon export than the other two models.The main conclusions of our study are that high-resolution interdisciplinary data from autonomous platforms have enormous potential for optimizing biological models, and that different models can fit surface observations almost equally well, yet yield very different estimates of permanent carbon export, emphasizing the importance of validating carbon cycle models in the twilight zone.

Fig. 1 .
Fig. 1.Location of the NAB08 experiment.Path of the Lagrangian float is shown in blue.R/V Knorr cruise track is shown in red.

Fig. 2 .
Fig. 2. Schematic of the 2p2s biological model.The 1p1s and 2p1s models are simplified variants of this model.
slope of the small phytoplankton P vs.I curve α Dia mmol N (mg Chl d W) −1 m 2 initial slope of the diatom P vs.I curve R Si:N mol Si (mol N) mmol N) −1 maximum Chl-to-N ratio of diatoms R Oxy:DIN mmol O (mmol N) −1 m is the number of observations per data type, W m is an arbitrary weight for each data type, x m,obs i is an individual observation of data type m, x m,mod i (p) is the corresponding model counterpart and depends on the vector of biological parameter values p.The values of N m and W m are given in Table 4.The weights can be interpreted as inverse variances for each data type multiplied by an arbitrary scalar and thus carry units of inverse variance.They were chosen to give roughly equal weight to the five data types.The assimilated variable types are Chl, PON, DIN, Oxy and Si, where PON = Phy + Dia + Det N + Zoo or, in case of the 2p2s model, Phy + Dia + Det N + Zoo + Cys.

Fig. 4 .
Fig. 4. Observed and simulated mean chlorophyll concentrations for different depth intervals.

Fig. 5 .
Fig. 5. Observed and simulated mean PON concentrations for different depth intervals.

Fig. 6 .
Fig. 6.Observed and simulated mean DIN concentrations for different depth intervals.

Fig. 7 .
Fig. 7. Observed and simulated mean oxygen concentrations for different depth intervals.

Fig. 8 .
Fig. 8. Observed and simulated mean silicate concentrations for different depth intervals.
of 22.5-50 mmol C m −2 d −1 and by Bury et al. (2001) of 40-47.5 mmol C m −2 d −1 for the upper 350 m, but much higher than the estimate by Lochte et al. (1993) of 9.8 mmol C m −2 d −1 at 150 m depth.

Table 1 .
Data used in modeling.All data are from the Lagrangian float except for silicate, which is from bottle samples.Each day the float profiled to about 250 m.For the rest of the day, it sampled continuously within the mixed layer.Estimated errors in a single 1-m by 1-day bin are shown for variables used in the optimization.
2 During daily profiles.

Table 2 .
State variables and parameters of the biological model.

Table 3 .
Range from which initial model parameter guesses are drawn, optimized parameters and their corresponding a posteriori errors as variance (from Hessian matrix analysis) scaled by parameter values.Errors greater than 100 % are highlighted in dark gray, errors greater than 50 % are highlighted in light gray.

Table 4 .
Number of assimilated observations per data type (N m ), cost function weights (W m ), optimal cost function values (totals and individual contributions from the different variable types) and associated standard errors resulting from uncertainties in the observations.The standard errors were calculated by Monte Carlo analysis as described in the text.

Table 5 .
Simulated POC export at 100 and 600 m (total export integrated over the simulation period and average flux) and contributions of each variable.
Data collected by other platforms, such as the seagliders and ship, Bagniewski et al.: Optimizing models of the North Atlantic spring bloom are not easily interpreted as Lagrangian timeseries and are thus difficult to use for validation.