Can biogeochemical fluxes be recovered from nitrate and chlorophyll data? A case study assimilating data in the Northwestern Mediterranean Sea at the JGOFS-DYFAMED station

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

One of the principal objectives of studying biogeochemical cycles is to obtain precise estimates of the main fluxes, such as total, new and export oceanic productions. Since models can incorporate the a priori knowledge of the most important processes, they are increasingly used for this purpose. However, biogeochemical models are characterized by a large number of poorly known parameters. Moreover, the available data are rather sparse in both time and space, and represent concentrations, not fluxes. Therefore, the major challenge is to constrain the relevant fluxes using information from a limited number of observations and from models incorporating poorly known internal parameters.
The present study attempts to meet this challenge. In a 1D framework at the DYFAMED station (NW Mediterranean Sea), near-monthly nitrate and chlorophyll profiles and daily surface chlorophyll concentrations are assimilated in a coupled dynamical -biological model using the tangent linear and adjoint models. Following sensitivity analyses that show that some parameters cannot be recovered from the data set used, assimilation of observed 1997 data is performed. The first inversion considered clearly shows that, in agreement with previous studies, (1) the data impose a C/Chl ratio that varies with depth (i.e. light) and (2) the ''initial'' conditions (e.g. winter nitrate profile) strongly constrain the annual biogeochemical fluxes. After assimilation of the 1997 data, the agreement between the data and the model is quantitatively improved in 1995 and 1996, which can be considered a good validation of the methodology. However, the order of magnitude of the biogeochemical fluxes, and especially of the particulate export and regenerated production, are not correctly recovered. An analysis of the simulations shows that this result is associated with a strong decrease in zooplankton concentrations. An additional constraint of maintaining acceptable levels of zooplankton is therefore added. The results are improved, but remain unsatisfactory. A final inversion, which takes into account the a priori estimates of the major annual fluxes, is then performed. This shows that there is no

Introduction
Together with ocean dynamics, the main biogeochemical fluxes such as total, new and export productions, control the total inorganic content of the mixed layer, which in turn affects carbon dioxide exchanges between the ocean and the atmosphere. The study of biogeochemical cycles endeavours to estimate these fluxes accurately. Modeling approaches have the advantage of being able to make use of a priori knowledge, and are increasingly used for such purposes. But the solutions of biogeochemical models depend to a large extent on the values of the biological parameters. These parameters are not only numerous, but most of them are also poorly known since they often represent unmeasurable variables used in the parameterization of complex exchanges between biological compartments. Moreover, these compartments represent the averaged behavior of different biological species. Traditionally, the parameters are tuned until the best agreement between model results and data is reached. Data assimilation techniques provide a powerful means to achieve this tuning, and are becoming a mandatory step in model development. However although theoretically appealing, the practical use of data assimilation is a difficult task, and exploratory studies are needed before data assimilation can be routinely used in 3D biogeochemical models.
The available biogeochemical data are rather sparse in time and space, and often represent concentrations and not fluxes. Indeed, biogeochemical fluxes are difficult to measure on a routine basis. Primary production measurements require in situ incubations, while export measurements are generally evaluated from the deployment of sediment traps. Moreover, the error associated with these measurements is quite high (Richardson, 1993;Buesseler et al., 1992Buesseler et al., , 2000. Nitrate and chlorophyll data, on the other hand, are more accessible, more widespread, and more precise.
The major challenge is therefore to try to constrain the biogeochemical fluxes using information from a limited number of stock observations and nonlinear models associated with poorly known internal parameters. Such an attempt is presented in this study. In a 1D framework at the DYFAMED station, near-monthly nitrate and chlorophyll profiles and daily surface chlorophyll concentrations are assimilated in a coupled dynamical -biological model using the tangent linear and adjoint models.
They are several reasons why the DYFAMED station, located in the Northwestern Mediterranean Sea, is an interesting test case. First, several biogeochemical production regimes that take place in the world ocean are found here. Deep convection occurs during winter, leading to a spring bloom. Oligotrophy prevails during summer while perturbations in the meteorological forcing generate a secondary bloom in fall. Secondly, the station is far enough away from the Ligurian Current to be sufficiently protected from lateral transport, thereby permitting a 1D study. Moreover, DYFAMED is a JGOFS time-series station, which means that a data set of biogeochemical and physical parameters is available to carry out and validate simulations. A final reason for using the DYFAMED station is that it is relatively well known (Deep-Sea Res., Part 2, special issue, 29 (12), 2002, in press), and has been the subject of previous model studies (Lévy et al., 1998a;Mémery et al., 2002).
Several inverse methods have been applied to the estimation of biogeochemical parameters. These methods are all based on the minimization of a cost function that depends on the parameters to be estimated and measures the distance between the observations and the corresponding values calculated by the model. But as pointed out by Athias et al. (2000), the result of the assimilation may depend on the choice of the inverse method used to find the minima of the cost function. It is thus of primary importance to choose an adequate method and test its performance. Most previous studies have used a variational data assimilation approach, based on a linearization of the cost function around an a priori solution (Prunet et al., 1996a,b), and a gradient descent method (Fennel et al., 2001;Gunson et al., 1999;Lawson et al., 1995Lawson et al., , 1996Schartau et al., 2001;Spitz et al., 1998Spitz et al., , 2001. Other studies have used nonlinear optimization techniques (Evans, 1999;Fasham and Evans, 1995) or global optimization algorithms such as simulated annealing (Athias et al., 2000;Armstrong, 1996, 1999;Matear, 1995).
Mathematical arguments favor the choice of nonlinear techniques. Indeed, the inverse parameter estimation problem is always a nonlinear problem, even in the case of linear dynamics (Evensen et al., 1998). Therefore, the cost function may possess different local minima. There are, however, three reasons why, in our opinion, the optimal control method, often referred to as the adjoint method, is promising. First of all, global optimization algorithms require a large number of cost function evaluations and this does not seem possible using a full three-dimensional model, which is the final goal of ocean color data assimilation. Secondly, the adjoint method is already being used with some success in ocean global circulation models, which is a first step towards its implementation in 3D biogeochemical models. Thirdly, the a priori parameter values or ''first guess'' values required in the adjoint method are not completely unknown, even if they are attributed large error bars. It therefore seems natural to treat the inverse problem locally around these a priori values. The ideas and mathematical concepts of optimal control theory were formalized about 30 years ago (Lions, 1971) and have since received a lot of attention with respect to their applications in meteorology (Le Dimet and Talagrand, 1986;Courtier and Talagrand, 1987;Talagrand and Courtier, 1987) and oceanography (Long and Thacker, 1989a,b;Luong et al., 1998).
For simplicity, most previous biogeochemical data assimilation studies have considered 0D model set-ups, in which biogeochemical variables corre-spond to mixed-layer values. This approach is suitable when most biogeochemical activity occurs within the mixed layer. This is the case in the spring phytoplankton bloom, provided that the upload of nutrients is correct. It is, however, not the case for oligotrophic regimes, where the chlorophyll maximum is below the mixed layer. The novelty of our approach is in the use of the adjoint technique in a 1D context to assimilate observed data, as well as the first use of the DYFAMED data set for data assimilation purposes.
The biogeochemical model used in this study comprises six nitrogen compartments. In a study where data from the Bermuda Atlantica Time Series (BATS) station are assimilated, Spitz et al. (2001) warn that data cannot be satisfactorily assimilated if the structure of the biogeochemical model is not adequate. We follow the procedure adopted by Spitz et al. (2001) and use the results of a first assimilation attempt to modify the model. We show that a significant improvement in model performance is obtained when the C/Chl ratio is allowed to vary in the vertical, and also when the winter nitrate upload is estimated.
The paper is organized as follows: after a presentation of the model and the method, sensitivity analyses are carried out to determine which parameters can be recovered from the data set used. Then an assimilation of observed data for the year 1997 is performed. After the assimilation, some of the fluxes (e.g. particulate export production) are shown to be very poorly estimated.
An analysis of the simulation shows that zooplankton values have fallen to a very low level, which would explain this result. An additional constraint is then added: the zooplankton stock during the post bloom regime (Andersen and Prieur, 2000). Although the results are improved, they are still not satisfactory. The method is then validated using independent data from years 1995 to 1996. A final inversion is then carried out that also takes into account the a priori estimates of the major annual fluxes (total and new productions, and particulate and dissolved exports). This inversion shows that there is no inconsistency between the nitrate and chlorophyll data, the estimates of the fluxes and the model used. The difficulty of recovering the biogeochemical fluxes from real data assimilation of concentrations and stocks (nitrate, chlorophyll, and eventually zooplankton) is then discussed.

Model description
The biogeochemical model used in this study is a NNPZD-DOM model (Nitrate, Ammonium, Phytoplankton, Zooplankton Detritus, Dissolved Organic Matter), qualitatively calibrated for DYFAMED by Mémery et al. (2002). Compared to other more complete models (Fasham et al., 1990;Lévy et al., 1998a), its relative simplicity is a trade-off to obtain a first approximation of the basic biogeochemical fluxes with a minimum number of prognostic variables. Nitrate and ammonium allow the estimation of new and regenerated production, zooplankton mortality and detrital sedimentation feed the particle export flux, and winter mixing of accumulated semi-refractory DOM is associated with the dissolved export flux, which can be the major export flux in the Mediterranean Sea (Copin- Montégut and Avril, 1993). A schematic representation of the model is shown in Fig. 1, and the basic equations are presented below. More details about the modeled biogeochemical processes can be found in Mémery et al. (2002) and in Lévy et al. (1998a).
The biogeochemical model is embedded in a 1D physical model, which simulates the evolution over time of velocity, temperature, salinity and turbulent kinetic energy (TKE). As for dynamical processes, the only one taken into account is vertical diffusion. The mixing coefficient, K T , is obtained diagnostically from TKE, with a 1.5 closure scheme in the Mellor Yamada nomenclature (Gaspar et al., 1990).T h e model covers the first 400 m of the water column, with a vertical discretization of 5 m.
The biogeochemical tracers are vertically mixed with the same diffusion coefficient as temperature and salinity. A specific reaction term, F C , is added to the diffusion equation. Tracers are expressed in terms of their nitrogen contents (mmol N m À 3 ). For each of the state variables, NO 3 ,N H 4 , P, Z, D and DOM, the prognostic equation reads as follows: where C is the tracer concentration. For the surface layers (0 -150 m), the reaction terms read as follows: The parameters are presented in Table 1. The formulation of phytoplankton growth takes into account limitation by both nutrients and light. Nutrient limitation follows the Hurtt and Armstrong (1996) kinetics: Based on the hypothesis that the total limitation follows the same law, i.e.
we obtain, Light limitation is expressed as follows: The photosynthetic available radiation (PAR) is predicted from surface irradiance and phytoplankton pigment content according to a light absorption model. Two different wavelengths are considered and the absorption coefficients depend on the local phytoplankton concentrations: PARðz; PÞ¼PAR r ðz; PÞþPAR g ðz; PÞð 12Þ PAR r ðz; PÞ¼PAR r ðz À Dz; PÞð1 À expðÀk r DzÞÞ where k go = 0.0232 m À 1 and k ro = 0.225 m À 1 are the water absorptions in green and red. Grazing of phytoplankton and detritus is formulated following Fasham et al. (1990): Other modeled biogeochemical interactions include phytoplankton mortality, phytoplankton exudation, zooplankton mortality (considered as large particles which are supposed to be instantaneously exported below the productive layer and remineralized in the water column), zooplankton exudation, fecal pellet production, detritus sedimentation, detritus breakdown, nitrification, and dissolved organic matter remineralization (Fig. 1).
Below a depth of 150 m (referred to as z bio hereafter), remineralization processes are preponderant and the surface model does not apply. Instead, decay of phytoplankton, zooplankton and detritus in nutrients, and a vertical redistribution of zooplankton mortality according to Martin and Fitzwater's (1992) profile parameterize remineralization below the surface layer. This parameterization conserves total nitrogen.
The remineralization of the export flux due to zooplankton mortality (large particles) is expressed as: Stating that no deposition occurs on the ocean floor, a further condition is added:

The DYFAMED simulation set-up
The model set up is the same as in Mémery et al. (2002). The standard run consists of a spin-up period of two repeated years (1995), followed by two more years, 1995 and 1996. The end of year 1996 is the restart point for the working year 1997. A complete description of seasonal and interannual variability is documented in Mémery et al. (2002). The simulation is forced with ECMWF atmospheric data, which give the wind stresses and heat fluxes every 6 h. Although the DYFAMED station can be considered, at a first approximation to behave like a 1D vertical system (Lévy et al., 1998a;Andersen and Prieur, 2000), deep convection in winter involves barotropic -baroclinic instabilities, which are three-dimensional processes (Madec et al., 1991;Lévy et al., 1998b). In winter, therefore these 3D processes must be parameterized. We use a crude parameterization consisting of restoring temperature, salinity and dissolved tracers (nitrate and DOM) towards homogeneous profiles in February when convection is at its peak and the mixed-layer is deepest. For example, the NO 3 reaction term, F NO 3 , becomes F NO 3 + a(N r -NO 3 ). N r is the restoring profile, and a equals 1/7 day À 1 in February and 0 otherwise. For temperature and salinity, we use the February data as restoring profiles, since they are fairly homogeneous. For nitrate, the only available winter profile is not homogeneous. Therefore, we construct a homogeneous profile with the a priori averaged nitrate value below 200 m (6 mmol N m À 3 ). This restoring profile can be seen as the nitrate initial condition or preconditioning, since it determines to a large extent the amount of nitrate made available to the system. For DOM, restoring is carried out towards zero, which implicitly assumes that all the semi-labile DOM is exported in winter through convection.
In this work, as well as in the previous studies carried out at DYFAMED (Lévy et al., 1998a;Mémery et al., 2002), advection is neglected. This might result in a crude approximation in summer during strong wind events. Indeed, Andersen and Prieur (2000) and Chifflet et al. (2001) suggest that Ekman pumping may be responsible for an upward advective shift of the nitracline and the deep chlorophyll maximum (DCM) of the order of 10 to 15 m.

The DYFAMED data set
The data used in this study are monthly chlorophyll and nitrate profiles collected for the years 1995 -1997 at the DYFAMED station, and surface fluorescence data measured by the Carioca buoy (Hood andMerlivat, 2001) moored at DYFAMED in 1997 (Figs. 2 -4).We also make use of various annual estimates of new production, regenerated production, dissolved and particulate export (estimated from disparate measurements undertaken in the 1990s), and of the zooplankton content estimated during a special cruise in May 1995 (Andersen and Prieur, 2000).
The evolution over time of the chlorophyll and nitrate profiles (Figs. 3 and 4) reflects the seasonal variability at DYFAMED. Winter mixing brings nutrients to the surface, but the short residence time  of algae in the euphotic layer, swept along by strong vertical motions, prevents the development of biomass. As the year progresses, the surface layer becomes more stable, thus allowing the winter nutrient enrichment to be utilized continuously. As a consequence, the algae bloom. From mid-May to November, the situation remains fairly stable with the upper layer nutrient content very low, and the system mainly oligotrophic and characterized by a deep chlorophyll maximum. It may be noticed that the nitrate profiles show strong variability below 100 m. This variability cannot be attributed to biological processes since these occur closer to the surface, and it is therefore most likely due to advection. It cannot therefore be captured by the model. During the oligotrophic period, the location of the nitracline and the DCM is fairly constant. This may be an indication of the absence of strong Ekman pumping and a crude justification for the neglect of vertical advection.
Estimates of total production (TP) vary considerably. Minas (1970), using in situ measurements, estimated the Northwestern Mediterranean TP to be about 78 g C m À 2 year À 1 , whereas Antoine et al. (1995) obtained a value of 157.5 g C m À 2 year À 1 from satellite data and Marty and Chiavérini (2002) an average value of 156 g C m À 2 year À 1 over the 1990s from C14 incubations (although they warn that this value should be considered as in the upper range). Estimates of new production (NP) range from 20 (Bethoux, 1989) to 42 g C m À 2 year À 1 (Marty and Chiavérini, 2002). This latter value is based on nitrate uptake (NU) and does not therefore correct for nitrification, which explains the difference with the Bethoux (1989) estimate. With regard to export fluxes, the DOC export flux has been estimated by Copin-Montegut and Avril (1993), based on 1991 DOC data, to be around 14.8 g C m À 2 year À 1 , while sediment trap data from 1987 (Miquel et al., 1993) predict a particulate export at 200 m of 4 g C m À 2 year À 1 .
The data assimilation period is restricted to 1997, when all three types of synoptic data are available.
For coherence with the model set-up strategy, the non-homogeneous February nitrate profile is not accounted for in the assimilation. Data from years 1995 and 1996 are kept for validation.
In order to compare model results and data, calculated phytoplankton concentrations are converted to mg Chla m À 3 using a constant (in time, space and during parameter optimization) Red field ratio C/N of r d = 6.625, and a C/Chl ratio (r c in Table 1) that will be optimized. A linear relationship is assumed between fluorescence and chlorophyll concentrations. The proportionality coefficient was chosen so that surface fluorescence data and monthly chlorophyll profiles were coherent with one another.

An initial simulation
A first run is carried out using the a priori parameter set and without any data assimilation . The seasonal variability is well reproduced, but the spring bloom appears to be too strong and the transition towards oligotrophy poorly represented. On the other hand, and because the parameters have been tuned for this purpose (Lévy et al., 1998a), this ''first guess'' run predicts reasonably good fluxes: 110.6 g C Fig. 5. Three smallest eigenvalues and corresponding eigenvectors of the Hessian for the 25-parameter optimization case. m À 2 year À 1 for TP, 22.5 g C m À 2 year À 1 for NP, 13.2 g C m À 2 year À 1 for DOC export and 5.2 g C m À 2 year À 1 for particulate export.

Control variables and cost function
The variational data assimilation method mainly consists in finding an optimal control minimizing a cost function which measures the distance, in a weighted least-squares sense, between the model's solution and the observations. In the current literature, authors usually consider the biological parameters as being the control variables. It is therefore assumed that the discrepancy between the model and the observations can be attributed mainly to these poorly known quantities. However, implicitly assuming that the model equations are perfect except for the unknown parameters may be a problem (Evensen et al.,1998). Firstly, minimization may lead to an optimal parameter set attempting to correct for errors other than those related to the particular parameterizations. Secondly, the tuning of the parameters may not be sufficient to fit real observations and it may be necessary to introduce more degrees of freedom in the model equations.
In our case, the winter restoring terms present in the equations are clearly crude parameterizations made to take into account 3Dprocesses which are not explicitly resolved. The restoring value of 6 mmol N m À 3 for NO 3 was chosen empirically on the basis of the DYFAMED data. The NO 3 restoring term is particularly important in the biological dynamics (Mémery et al., 2002). It can be regarded as an initial condition from which we want to start at the end of winter. It can also be regarded as a particular form of a model error term, and we may need to control it.
The 1D formulation of the model as well as a first real data assimilation experiment (Case 1, Section 4.4.1) also lead us to consider a C/Chl ratio r c that varies with depth. We have therefore considered different cases in this study. Case 1 is classical: the cost function depends on the biological parameters, which are undistributed, or constant in time and space. In Case 2, the cost function depends on all the previous parameters but r c can vary with depth. It also depends on the NO 3 restoring profile, N r , which is also space-distributed.
As shown in Table 1, the biological parameters have very different orders of magnitude. To avoid any numerical difficulties which might arise from this during the minimization, we adimensionalize the parameter vector K, dividing each parameter K i by its first guess value K i 0 , In the same way, we consider the dimensionless NO 3 restoring profile, n r =(n i r ) i =1,...nz = N r /6, where nz is the number of space discretization points. Such a non-dimensionalization procedure can be regarded as a preconditioning for minimization.
To sum up, in Case 1, the control variable is x = k of size p, where p is the number of parameters, and k is dimensionless. In Case 2, the control variable is x=(k,r c ,n r ) of size ( p À 1) + nz + nz, where k represents the parameter vector without r c and n r . In both cases, the first guess value for x is x 0 =(1,...,1).
The model-data misfit part, J 0 , of the cost function can be written as the sum of three terms:  We also add penalty terms to J 0 . The first term, accounts for the a priori parameter values and their standard deviations, r i ¼ 1 ffiffiffi ffi w i p , given in Table 2. This term enables the minimization to avoid biologically absurd optimal parameter values.
The second term, which only applies in Case 2, is a Tikhonov regularization term (Tikhonov and Arsenin, 1977).
This is added to avoid instability phenomena inherent in such inverse problems. Its role is to smooth the space-distributed control n r . The choice of the weight w is not straightforward, as a compromise has to be found between numerical stability and the quality of the minimization of J 0 . In practice, w was chosen so that the data misfit part of the cost and its regularization part are balanced. The same regularization term is added for the distributed parameter r c .

Automatic differentiation
To minimize the cost function J, we applied the quasi-Newton algorithm implemented in the n1qn3 Fortran subroutine of Gilbert and Lemaréchal (1989). The computation of the gradient of J with respect to control variables is required at each step of the minimization. This gradient results in one integration of the adjoint model. The tangent and adjoint codes were obtained using the automatic differentiation program Odyssée (Faure and Papegay, 1997;Griewank, 2000), which is an efficient tool for obtaining adjoint codes since it enables the automatic production of adjoint instructions. However, codes produced by automatic differentiation do not usually use computer memory in a very efficient way. Saving the direct model trajectory is the major problem. A differentiation program has to follow systematic methods to provide the evaluation trajectory. Thus, Odyssée systematically uses a local calculation and storage technique for the trajectory. Automatically differentiating a 3D model and using the adjoint code directly seems impossible. This problem partially disappears using a 1D model, so that the work to be done on the automatically produced adjoint routines is more reasonable. In this study, we used an integration period of 1 year, the choice being made to cover a whole annual cycle of biological activity when all three types of data were available. Moreover, it involves affordable computer memory resources and acceptable calculation times. An adjoint run is about 10 times more expensive than a direct run in terms of memory requirements and CPU time.

Sensitivity analysis method
As was highlighted in Fennel et al. (2001),a sensitivity analysis should be an integral part of any attempt to optimize the parameters of a biological model. Fennel et al. identified different inconsisten- cies in the formulation of their inverse problem such as parameters which could not be determined independently.
The a priori approach consists in using a sensitivity analysis based on twin experiments, that is in the case where the global minimum of the cost function is known, to help determine the most sensitive and consequently the most important parameters to optimize in order to solve a numerically well-posed inverse problem.
Let us consider that the cost function only depends on parameter k and that no penalty term is added: A small perturbation h on parameter set k induces a perturbation on /, /ðk þ hÞÀ/ðkÞc/VðkÞh /V(k) is the sensitivity matrix, of size N obs Â p.I t can easily be calculated by p integrations of the tangent linear model, perturbing the parameters one after another. The ith integration leads to the ith column of matrix /V(k) that represents the influence of parameter k i on the model's output used to calculate the cost function. Then, if we expand J in its Taylor series about the optimal solution k 0 and ignore terms of a higher order than two, we obtain the following, as the first derivative vanishes: In the context of twin experiments, where data are generated by a model run using k 0 as the parameter set, the formula reduces to where H = /V(k 0 ) T W/V(k 0 ) is the Hessian of the cost function in this particular case. H is easily calculated once /V(k 0 ) has been obtained. It may be noted that H does not depend on the quality of the data but only on the intrinsic properties of the model structure, or in other words on the way in which the parameters are mutually related in the conceptual representation of the biological phenomena. The Hessian matrix provides a good approximation of the covariance matrix for the model parameters (Thacker, 1989). Its diagonal element h ii gives the sensitivity of the cost function to change in parameter k i , whereas its off-diagonal elements indicate the degree to which pairs of model parameters are correlated. In fact, all the needed information is contained in the eigenvalues and vectors of H. The condition number, k max /k min gives an indication of the degree of singularity in the problem. Small eigenvalues correspond to large uncertainties in the identification of the parameters that make a significant contribution in the related eigenvectors. Such an analysis provides a systematic method for determining parameters which may be difficult to identify.

Description of the experiments
This short subsection briefly presents all the different numerical experiments which follow in the text. Let us first remember that in Case 1 the cost function is assumed to depend on constant parameters, whereas in Case 2 it is assumed to depend on constant parameters and on two parameters (n r and r c ) allowing vertical variability. We first conduct a sensitivity analysis based on a twin experiment identical to Case 1 (Section 4.2). A synthetic data set is created with the model, by running the reference 1997 simulation with the parameter set k 0 . These synthetic data have the same spatio-temporal distribution as the real 1997 data. We then conduct twin experiments for both Case 1 (Section 4.3.1) and Case 2 (Section 4.3.2). Once these validation steps have been successfully completed, we come to the assimilation of real (i.e. observed) data. Section 4.4.1 deals with Case 1 optimization. Case 2 is divided into Case 2a in Section 4.4.2 and Case 2b in Section 4.4.3. As explained in Introduction, in Case 2b zooplankton data is added to the cost function, and in Case 2c annual flux estimates are also assimilated.

Twin experiment sensitivity analysis
We apply the sensitivity analysis methodology described in Section 3.3 to the twin experiments to ensure that we are not trying to solve an ill-posed inverse problem. Within this framework, the cost function is defined without any penalty term (Eq. (28)). The Hessian, H, is calculated at the optimal point k 0 , and then its eigenvalues and eigenvectors.
If the cost function is assumed to depend on the 25 parameters in Table 1, the condition number of its Hessian is 1.149 Â 10 17 , indicating that the inverse problem is extremely ill-conditioned. The three smallest eigenvalues, k 25 = 8.804 Â 10 À 15 , k 24 = 5.405 Â 10 À 5 and k 23 = 2.061 Â 10 À 4 and their corresponding eigenvectors, v 25 , v 24 and v 23 are plotted in Fig. 5. v 25 has contributions mainly for k rp , k gp and r pg . As the associated eigenvalue is almost 0, these parameters cannot be determined with the data under consideration. All of them belong to the optical model formulation (Eqs. (12) -(16)). Eigenvalue k 24 = 5.405 Â 10 À 5 is still quite small, and the important contributions in v 24 again correspond to the optical parameters. These parameters enter the model only in the combinations k gp (12r d / r pg r c ) lg and k rp (12r d /r pg r c ) lr (Eqs. (15) and (16)). It is therefore not surprising that these parameters appear to be dependent on one another. As a considerable improvement is needed in the condition number of H, we will not try to estimate the five optical parameters in all of the following experiments. Their values are kept fixed at the initial guesses. Only the C/Chl ratio r c will therefore play a role in the calculation of light limitation L I during the optimization procedure.
Under this new assumption, the condition number of the Hessian matrix falls to 7.489 Â 10 6 which is much more acceptable and comparable to the condition number values presented in F e n n e le ta l . (2001). The worse resolutions are observed for h r , c and m z . The coefficient for Martin's remineralization profile, h r , only concerns the remineralization model (below 150 m), and is kept constant in the following sections.
Finally, the chosen formulation includes 19 parameters. The condition number of the Hessian is 9.108 Â 10 5 . Eliminating two or three additional parameters, such as c, m z or k n , which are the least sensitive (Fig. 6), does not considerably improve the conditioning of the problem and may, on the con- trary, deprive the model of useful degrees of freedom.

Case 1
The second and final validation step before assimilation of observed data is to conduct twin experiments. Synthetic data are produced by the model using the first guess parameter vector k 0 containing the 19 selected parameters. In order to fully test the possibility of recovering these parameters from the synthetic data, no penalty term is added. The variances assumed are as follows:r cs = 0.3 mg Chla m À 3 , r cp = 0.1 mg Chla m À 3 , and r np = 0.4 mmol N m À 3 . This provides a good balance between the three terms of J 0 . Different first guesses for the parameter vector were obtained by perturbing k 0 within the a priori error deviation ranges given in Table 2. All the corresponding optimizations converged to the minimum of the cost function, J(k 0 )=0.
The results of this experiment are shown in Fig. 7. Table 2 presents the first guess parameter set k 0 + dk. The convergence criterion (NjJN V e, where e is a small value) is satisfied after 129 iterations. A cost function value of 7.35e À 11 is obtained, indicating that it has reached its global minimum. All the 19 parameters have been recovered. In the first 20 iterations, the cost function decreases rapidly then enters a long period of more gradual decline. This corresponds to the rapid recovery of the most sensitive parameters, followed by a difficult search for the optimal values of the less sensitive parameters (k n , c, and m z ).

Case 2
We are now concerned with the optimization of 18 space-and time-constant parameters and 2 space-dependent parameters, r c (the C/Chl ratio) and n r (the NO 3 winter preconditioning). This greatly increases the number of control variables from 19 to 178. However, the problem still remains overdetermined as the control dimension does not exceed the number of data. Nevertheless, this does not mean that it is possible to recover the model controls given the set of data we have for year 1997 at the DYFAMED station. Some controls may not have a real influence everywhere on the simulation results. Moreover, the method has to be able to distinguish between wrong biological parameter values, wrong C/Chl ratios and wrong NO 3 preconditioning to explain the discrepancy between Fig. 7. Convergence of the 19 parameter set from k 0 + dk to k 0 . model and data. To check this, a twin experiment is conducted in which no penalty term is added to the cost function relating to the constant parameters. A Tikhonov regularization term is used for the two distributed controls r c and n r .
J ðxÞ¼J 0 ðxÞþJ r ðn r ÞþJ r ðr c Þ Exact equivalents of year 1997 DYFAMED data were created using the reference a priori k 0 vector for the 18 parameters, the distributed r c ratio and the distributed n r profile plotted in Fig. 8. First guess control variables consist of a profile n r homogeneous on the vertical with value n r = 1 and the parameter set k 0 + dk in Table 2, except for r c which is also chosen to be homogeneous in the vertical. The first guess for r c is 1. Below 200 m, this value is the same as that used to generate the data. This choice was made because below 200 m r c can take any value without influencing the simulation or, therefore, the cost function value. The twin experiment proved to be successful, with all biological parameters, as well as r c , completely recovered and the restoring profile, n r , also very well reproduced. The r c and n r used to generate the data, first guess and recovered values, are plotted in Fig. 8. It takes 543 iterations for the optimization to converge. This increase in the number of iterations required compared with Case 1 is due to the increase in the size of the control vector.

Estimation of constant parameters-Case 1
A real data assimilation experiment is performed using nitrate and chlorophyll data from the DYFAMED station. The first guess is chosen as k 0 , which corresponds to the parameter values Fig. 8. Left: r c , first guess (long dashes), reference (continuous line) and optimal (short dashes superimposed on the reference). Right: n r , first guess (long dashes), reference (continuous line) and optimal (short dashes). The reference values, R c = 55 mg C mg Chla À 1 and N r = 6 mmol N m À 3 , are used for adimensionalization. commonly used in this area (Lévy et al., 1998a;Mémery et al., 2002), and the data used are those described in Section 2.2. A penalty term is added to the cost function of the twin experiments: This regularizes the inverse problem and the optimization converges in 33 iterations. The optimal parameter set, k opt , is given in Table 3, and the overall quality of the optimization is quantified by the normalized cost related to each type of data in Table 4.
Correlation coefficients have also been computed to measure the shift between the model and the data. For surface chlorophyll, a time correlation coefficient is defined by: where / opt cs and dcs denote the mean values. In the same way, space correlation coefficients are computed for each profile.
The largest change in the parameter set is obtained for r c , the C/Chl ratio, which is used to convert phytoplankton into chlorophyll units. While the first guess for r c is 55 mg C mg Chl À 1 , the optimal value i s1 3 5m gCm gC h l À 1 .Other important changes concern parameters related to detritus, and thus with the export. Detritus breakdown rate, l d , decreases from 1.04e À 6 s À 1 to 4.524e À 7 s À 1 , while detritus sedimentation speed, v d , increases from 5 to 8.7 m s À 1 . These changes in parameters have the effect of  decreasing the chlorophyll concentrations and of increasing the export. Globally, the assimilation resulted in a considerable improvement in the model/data match, compared to the first guess experiment (Table 4). Bloom intensity is divided by a factor of three, in agreement with the observations (Fig. 2), and the time correlation coefficient increases from 0.87 to 0.92, indicating that the starting date and bloom duration are simulated more accurately. Chlorophyll profiles during the bloom (dates 02/23, 03/02 and 03/21 in Fig. 3) also show a better fit with the data. These improvements are mainly associated with the change in the C/Chl ratio. As with chlorophyll, nitrate is also better simulated and this is particularly apparent above the nitracline (Fig. 4). Nitrate concentrations at the surface were overestimated, and are now reduced. The reduction in nitrate concentrations is essentially achieved through the change in parameters l d and v d , and the associated increase in export. It may be noticed that the cost with respect to nitrate is globally higher than that relative to chlorophyll (Table 4). This is due to the poor fit between simulated and observed nitrate concentrations below the nitracline, where advective events seem to be important.
Despite the improvements noted above, however, certain problems remain. First, there is no improvement in summer, during the oligotrophic regime. The model has difficulties in reproducing the subsurface chlorophyll maximum. On 09/01, the space correlation coefficient falls from 0.865 to 0.061, indicating a substantial shift between observations and model results in terms of subsurface maximum depth. Second, the exported flux is too large by almost one order of magnitude (Table 5).
To address these shortcomings, the model is modified as follows. In order to better represent the chlorophyll subsurface maxima, the C/Chl ratio is allowed to vary on the vertical. Indeed, this ratio is known to be highly dependent on the photosynthetic available radiation (PAR) (Geider et al., 1997). In the Mediterranean Sea, it can vary from 30 to 200 mg C mg Chl À 1 . Moreover, earlier optimization studies have already emphasized the need for a variable C/ Chl ratio Armstrong, 1996, 1999;Schartau et al., 2001;Spitz et al., 2001). The amount of NO 3 present in the model is very dependent on the winter preconditioning profile (Mémery et al., 2002).I n order to improve the match between observed and simulated NO 3 profiles without excessively increasing the export, a space-distributed NO 3 winter preconditioning profile is introduced in the control vector.

Estimation of parameters allowing vertical variability-Case 2a
Using the space-dependent r c and n r described above as control variables improves all three parts of the cost J 0 (Table 4). Compared to Case 1, the summer chlorophyll profiles are better represented and improvements are also noted during other periods of the year (Fig. 3). Nitrate profiles are also in better agreement with the observations (Fig. 4). Fig. 9 shows the retrieved r c and n r profiles. The r c ratio is high (up to 143 mg C mg Chl À 1 ) near the surface, and rapidly decreases with depth (45 mg C mg Chl À 1 at about 40 m). Due to the regularization term, J r (r c ) in the cost function, r c remains close to its first guess value (55 mg C mg Chl À 1 ) below 50 m. This decrease in r c with depth is realistic and physically related to the interdependence of r c and PAR. Obviously, allowing r c to vary only with depth is an oversimplification, since it does not account for the seasonal variations of PAR. However, this experiment clearly shows that the approximation already greatly improves the results in relation to the constant r c case.
As for the optimal nitrate preconditioning profile, this is approximately homogeneous in the first 200 m of the water column (Fig. 9), consistent with the fact that the winter mixed layer is almost homogeneous down to this depth. It contains less nitrate than the Table 5 Flux data (from Marty and Chiavérini, 2002)  Values are given in g C m À 2 year À 1 . NU stands for NO 3 uptake, NP for new production (NU corrected from nitrification), TP for total production, PF for total particulate export flux (SF, the sedimentation flux + ZF, the zooplankton mortality flux) and DOMF for DOM export flux.
first guess (4.9 compared to 6 mmol N m À 3 ), which again indicates that there was too much nitrate in the first guess experiment. The optimal parameters are given in Table 3. They do not differ greatly from those of Case 1, but it should be noted that the parameters which were thought to vary because of the presence of too much NO 3 at the surface (l d , v d ) have optimal values closer to their first guess value than was the case in Case 1. In other words, the optimal values of these parameters do not suffer any more from the excess of nitrate, and we may expect them to be more representative of reality.
Regarding the fluxes, DOC export and NP are barely modified by the assimilation, and keep their acceptable first guess value. However, the detritus flux is still too high, and TP too low (Table 5).T P has fallen to 59.6 g C m À 2 year À 1 , which is well below the data range. The detritus flux is more than three times too high (12.6 g C m À 2 year À 1 ) and consists only of the sedimentation of slow sinking particles. Export through sedimentation of fast sinking particles (ZF), which is parameterized through the zooplankton mortality term, has almost vanished, whereas it represented one third of total particulate export in the first guess experiment.
The question arises as to why the model chooses to favor sedimentation of slow sinking particles over zooplankton mortality fluxes. Globally, in Cases 1 and 2a, zooplankton concentrations, which are not constrained by any data, are much weaker than in our first guess run (Table 6). This explains the fall in ZF, which behaves like Z 2 . Sensitivity studies show that the measurements, which cause the zooplankton concentrations to decrease, are the surface values of the chlorophyll profiles corresponding to the dates 03/21 and 04/15 (the two squares corresponding to days 80 and 105 in Fig. 2). These values impose a slow decrease in surface chlorophyll after the bloom maximum. This decrease is much more rapid in the first  19 guess run. The data assimilation forces the model to simulate relatively high surface chlorophyll values until the month of April. This is achieved through a strong decrease in grazing, mainly resulting from the decrease in zooplankton concentrations. In Case 1, the model runs right between the two abovementioned post-bloom data (Fig. 2) and ZF c 10 À 3 gCm À 2 year À 1 , whereas in Case 2a the model runs below these measurements and ZF c 10 À 2 gCm À 2 year À 1 . Given the strong influence of these two post-bloom measurements on the assimilation results, it may intuitively be thought that the sampling frequency at this period must play an important role. However, this cannot be easily tested and generalized in a rigorous manner. Twin experiments could be performed but this was decided against as one of the objectives of the study was to use the available DYFAMED data set.

Adding zooplankton information-Case 2b
The above observations all point to the need to constrain zooplankton concentrations in some way. Although no synoptic zooplankton data is available, the zooplankton content was estimated to be 20 mmol Nm À 2 in May 1995 (Andersen and Prieur, 2000). Simulated zooplankton concentrations in May in the first guess experiment are consistent with this value, but they are too low in Case 1 and Case 2a (Table 6). We therefore added this zooplankton data to the cost function and conducted a Case 2 type optimization, referred to below as Case 2b.
As expected, zooplankton concentrations in May (17.9 mmol N m À 2 ) are now more acceptable. The fit between observed and simulated chlorophyll and nitrate data is not altered much, except during the post-bloom period (day 04/15, Fig. 3). Indeed, as discussed above, the increase in zooplankton results in an increase in the grazing pressure on phytoplankton, particularly during the post-bloom. However, the results are still much better than in Case 1 ( Table 4).
The estimated C/Chl ratio, r c , and NO 3 restoring profile, n r , are very similar to those of Case 2a (Fig.  9). Observation of the other parameters shown in Table 3 also reveals that modifications are only slight compared with Case 2a. The most significant changes concern the detritus sedimentation speed, which decreases from 7.8 to 5.8 m day À 1 , the assimilated fraction of phytoplankton, which increases from 63% to 72%, and the nominal preference for phytoplank-ton, which decreases from 0.67 to 0.63. Zooplankton grazing rate increases by only 2%, while the zooplankton loss rates decrease by less than 1%. This behavior emphasizes the strong nonlinearites in the zooplankton equation. Improvement in the zooplankton content benefits the main fluxes (Table 5). The sedimentation of slow sinking particles (SF) decreases in response to the decrease in v d compared to Cases 1 and 2a. This leads to a more realistic value for total particulate export (EF), which now includes a weak contribution from fast sinking particles (ZF). TP has increased to an intermediate value between that of Case 2a and the first guess experiment.

Validation against independent data
To validate the methodology, the model is run for years 1995 and 1996 with the optimal parameter set and the C/Chl profile obtained after assimilation of year 1997 data (Case 2b).
The inter-annual variability in winter conditions is taken into account by not using the NO 3 preconditioning profile estimated for 1997, to simulate years 1995 and1996. Instead, for each of the years 1995 and 1996 we set the parameters to their values observed in Case 2b and only estimate the NO 3 preconditioning profile. This results in a significant improvement of the model fit to chlorophyll profiles for years 1995 and1996 (Table 7, Figs. 10 and11). Surface and subsurface maxima in spring and summer are globally better simulated, although several profiles in summer -fall 1995 (09/14/95, 10/05/95, 11/15/95) are not really  (Table 4), the agreement between data and model is globally better in 1996, but not as good in 1995. This might be linked to the ocean dynamics. Atmospheric forcing is characterized by much lower high frequency variability in 1996 than in the other years, which means that the evolution of the tracers over time is smoother (Mémery et al., 2002): vertical movements due to Ekman pumping or other dynamical processes, which are not considered in these simulations, could be weaker in 1996. The quantitative improvement of the simulation results in 1995 and 1996 after assimilation of the year 1997 can be considered a validation of the methodology.

Stock data and flux values-Case 2c
The series of numerical experiments conducted in this work shows that the principal biogeochemical fluxes can be significantly modified by the assimilation of stock data. In addition, the experiments illustrate the impossibility of recovering correct fluxes from the assimilation of NO 3 and chlorophyll data alone (Case 2a). When a single zooplankton measurement is added (Case 2b), flux prediction is improved. However, this improvement is not sufficient, and leads to a slight deterioration of the fit to chlorophyll data.
One possible reason for the lack of success in retrieving the fluxes could be that the model structure prohibits solutions consistent with both the stock data and the order of magnitude of the fluxes. To test this hypothesis, a final experiment was conducted in which all types of data were assimilated: stock measurements as well as the data on the four fluxes presented in Table  5, taken from Marty and Chiavérini (2002).T h i s experiment, referred to as Case 2c, proved to be successful. Surface chlorophyll, chlorophyll and NO 3 profile optimization results are comparable to those obtained earlier (Table 4), while the zooplankton con- Fig. 12. NO 3 profiles (mmol N m À 3 ) for year 1995. Depth in meters. The continuous line is the observation, the thin dashed line is the first guess run and the bold dashed line is the Case 2b optimal parameter set run. tent in May is 22.7 mmol N m À 2 , which is the best we obtain with the data (Table 6). More importantly, the predicted fluxes are now correct (Table 5). The hypothesis concerning the role of the model structure in the lack of success in retrieving the fluxes can therefore be rejected.
The estimated C/Chl ratio and NO 3 preconditioning profile are very similar to those found in Case 2a and Case 2b. With regard to the other parameters, presented in Table 3, the most striking changes compared with the previous cases appear for k n , c, l d and v d . The halfsaturation constant, k n , and the exudation fraction, c, were found to be among the less sensitive parameters in Section 4.2 and, not surprisingly, were not greatly modified by the data assimilation in Case 1, Case 2a and Case 2b. This was no longer true in Case 2c where production flux data seem to constrain these parameters. Sedimentation speed, v d , for which the optimal value was always found to be higher than the first guess value, is now lower, an observation that can be related to the decrease in the export flux imposed by the data. Furthermore, detritus feeds the NH 4 and DOM compartments at a higher rate l d .
Unlike for chlorophyll data, the relative cost of the NO 3 profiles is higher in Case 2c than in Case 2b (Table 4). Unlike in the first guess, the flux that had to be modified in Case 2c was particle export production (Table 5). To achieve this new balance, the inversion degraded the fit of the NO 3 profiles. Although a bias due to dynamical processes cannot be entirely discounted, this finding could result from an underestimation of particle export production, obtained from shallow sediment traps (Buesseler, 1991).

Conclusions
The present study has demonstrated the feasibility of using the adjoint data assimilation method to constrain the biological part of a one-dimensional coupled physical -biogeochemical model with observed data from the DYFAMED station in the Northwestern Mediterranean Sea.
A preliminary sensitivity analysis showed that not all parameters can be estimated and enabled the inverse problem to be correctly formulated. The correctness of the formulation was demonstrated using identical twin experiments.
Real surface chlorophyll data as well as chlorophyll and NO 3 profile measurements were assimilated for the year 1997. The assimilation method required choosing a judicious first guess. Different assimilation exercises then led to a number of conclusions.
As with earlier modeling studies (Geider et al., 1997), it was quantitatively shown that it is impossible to simulate correct surface chlorophyll bloom intensity together with a correct summer subsurface chlorophyll maximum in an oligotrophic regime using a constant C/Chl ratio. We used a C/Chl ratio allowing vertical variability to improve data assimilation results.
It was also quantitatively demonstrated that in our 1DDYFAMED context, the optimization of winter NO 3 conditions is of prime importance to data assimilation efficiency. The winter restoring profile governs the quantity of nitrate present in the model and needs to be controlled. Some estimated parameter values appeared to be dependent on this quantity. Moreover, NO 3 concentrations below the surface layer are largely imposed by this restoring profile.
Finally, a good improvement in the chlorophyll model-data fit was obtained. As for NO 3 data, the improvement was less significant below a depth of 100 m where data are affected by 3D physics not taken into account in the model. The method was then validated using an independent data set for years 1995 and 1996, and a new simulation conducted. This calculated acceptable fluxes and also predicted chlorophyll and NO 3 profiles more correctly than the first guess simulation.
The last but perhaps most important point concerns the simulation of production and export fluxes. It was shown that, using a simple NNPZD-DOM model, it is not possible to recover correct fluxes simply through the assimilation of NO 3 and chlorophyll stock data alone, although the model was able to simulate correct stocks together with correct fluxes. This failure was due to the collapse of zooplankton in the model. Nonlinearities, particularly in the zooplankton conservation equation, are certainly responsible for this behavior: minor changes in the parameters can considerably modify system evolution. Nevertheless, the addition of zooplankton data to the cost function did not prove to be entirely successful or enable the correct fluxes to be recovered. As stocks result from a balance between input and output fluxes, it may be that information associated with stocks alone cannot constrain the different fluxes independently. This conclusion implies that it may be impossible to estimate the biogeochemical fluxes accurately if a first guess of these fluxes is not taken into account and assimilated. Further studies, using other models and different data coverage, are needed to determine whether this statement can be generalized.