Optimization and sensitivity study of a biogeochemistry ocean model using an implicit solver and in situ phosphate data

A new implicit method for obtaining equilibrium solutions and their sensitivity to changes in parameters is described and applied to an OCMIP‐2 type ocean‐biogeochemistry model. The method is used to optimize model parameters by minimizing the difference between the observed and simulated PO4 distribution. The optimized parameters include (1) the exponent α in the power law vertical profile for particulate organic matter (POM) fluxes, (2) the fraction σ of biological production allocated to dissolved organic matter (DOM) and (3) the rate constant κ for the remineralization of DOM. Global PO4 observations constrain σ and κ but not independently because their sensitivity patterns are highly correlated. In contrast, the sensitivity pattern for α is uncorrelated to those of the other parameters, allowing it to be independently constrained. We show that export production from POC is well constrained by the distribution of PO4 in an OCMIP‐2 type model, but that new production and export production from DOC are not well constrained. With the optimal parameter set (α = −1.0, σ = 0.74, and κ = 1.0 yrs−1) the fraction of the spatial PO4 variance captured by our model increases from 60% with the reference OCMIP‐2 parameters to 70%. Combined changes in σ and κ account for most of the improvements by reducing but not completely eliminating the nutrient trapping effect in the Eastern Equatorial Pacific and northern Indian Ocean that causes the model to over‐predict PO4 concentrations. Important remaining model‐data misfits in the deep North Atlantic where PO4 is over predicted and in the North Pacific where the model does not produce the observed sharp nutricline are likely attributable to deficiencies in ocean transport. The fact that the fraction of unexplained variance is large at the optimal parameter values highlights the importance of properly simulating physical transport for ocean biogeochemical modeling.


Introduction
[2] The ocean is the largest dynamic reservoir of carbon on timescales of months to millennia. The wide range of timescales with which the ocean can interact with the atmosphere is a consequence of ocean transport that allows some water parcels to be ventilated relatively quickly while others are transported far from the surface and can remain shielded from the atmosphere for several thousand years. The oceanic carbon reservoir's broad range of response times makes it a particularly challenging system to model. For example, simulating the global carbon cycle over the relatively short anthropocene requires that the model be first spun up for several thousand years in order to distinguish trends associ-ated with anthropogenic perturbations from the model's drift toward its preindustrial equilibrium. Furthermore, the lengthy spin-up integrations must be repeated each time a model parameter is changed even if the parameter is associated with the biogeochemical model and does not lead to any change in the ocean circulation. As a result there has been few systematic parameter sensitivity studies of global biogeochemistry models. Parameter uncertainty thus remains large and poorly known.
[3] To overcome the computational costs associated with model spin-up, we have implemented a new solver for obtaining equilibrium solutions to a global ocean-biogeochemistry model. The new offline model makes use of timeaveraged ocean flow and eddy-diffusion tensor fields obtained from a dynamical ocean general circulation model (OGCM) and computes the biogeochemistry model's equilibrium chemical tracer distributions using an iterative approach based on Newton's method that avoids explicitly time stepping the governing equations. With the new solver, global three-dimensional equilibrium tracer distributions are obtained in minutes on a single processor workstation whereas the same solution obtained using the traditional time stepping approach would take weeks. In the present study, we take advantage of the new efficiency to perform an extensive parameter optimization and sensitivity study of an ocean-biogeochemistry model based on the formulation used for phase 2 of the Ocean Carbon Model Intercomparison Project (OCMIP-2, R. Najjar and J. Orr, Design of OCMIP-2 simulations of chlorofluorocarbons, the solubility pump and common biogeochemistry, 1998, http://www.ipsl.jussieu.fr/ OCMIP/phase2/simulations/design.ps) hereinafter referred to as Najjar and Orr online document). Our goal is to determine optimal parameter values along with their associated uncertainty.
[4] Our study is complementary to the goals of OCMIP-2 in that we study the sensitivity of changes in the biogeochemistry model with a fixed ocean circulation, whereas OCMIP-2 [e.g., Sarmiento et al., 2000;Matsumoto et al., 2004;Doney et al., 2004] used a fixed and common biogeochemistry model to focus on the impacts of changes in the ocean circulation as manifested in the suite of participating ocean models. Our study thus builds on the previous works of Anderson and Sarmiento [1995] and Tajika [1996, 1997] that studied the impact of changing export production and remineralization profiles on various chemical tracers in models that were structurally similar to the OCMIP-2 formulation. Our contribution is to present a much more extensive and systematic parameter sensitivity study.
[5] In our model, as in all OCMIP-2 type models, the sink and source of dissolved inorganic carbon (DIC) due to the production and remineralization of organic matter is keyed to the phosphorus cycle. Model parameters associated with the cycling of phosphorus between phosphate (PO 4 ), dissolved organic phosphorus (DOP) and particulate organic phosphorus (POP) are: (1) s, the fraction of net primary production going into DOP, (2) k, the first-order decomposition rate constant of semilabile DOP and (3) a, the exponent of the power law for the remineralization of POP with depth. In the present study, we report results dealing with phosphate and DOP. Results concerning the distribution of DIC and the air-sea exchange of carbon-dioxide will be presented in a follow up paper.
[6] Previous research aimed at optimizing large-scale ocean biogeochemistry model parameters have focused on the adjoint method [Matear and Holloway, 1995;Schlitzer, 2000Schlitzer, , 2002Usbeck et al., 2003]. Although there is some overlap between the goals of these studies and ours there are some important differences and these are reflected in the different methodologies. The common element between the previous studies and ours is that we try to estimate uncertain parameters by minimizing a cost function. The inputs and importantly the number of inputs to the cost functions are very different however. In the adjoint studies the inputs to the cost function include among other parameters the three-dimensional phosphate distribution, the circulation of the ocean and either a few biogeochemical model parameters in the case of Matear and Holloway [1995] or a very large number of parameters including surface export production in the case of Schlitzer [2000Schlitzer [ , 2002 and Usbeck et al. [2003]. The end result of the adjoint studies is a weighted average of the observations and the simulated fields. Because the number of adjustable parameters in these studies far exceeds the number of degrees of freedom in the observational data set, the uncertainty of the results cannot easily be quantified [Schlitzer, 2000]. In contrast, the inputs to our cost function include only three parameters related to the OCMIP-2 biogeochemical model. Because our cost function has only three dimensions we are able to plot its shape to determine how well the data constrains the parameters and thus quantify parameter uncertainties.
[7] The difference in the number of model parameters also underlies the differences in computational approach. Studies that involve cost functions with a large number of adjustable parameters are often based on adjoint methods because it allows for the efficient computation of the gradient of the cost function. The number of inputs in our cost function is kept to a minimum by holding the circulation fixed and considering only phosphate distributions that are equilibrium solutions of the model. This is made computationally feasible by the use of Newton's method to efficiently compute equilibrium solutions.
[8] The paper is organized as follows: In section 2 we describe how we use Newton's method to find equilibrium solutions. In section 2.1 we describe the OGCM simulation used to obtain the velocity and diffusivity fields needed to construct the transport operator. In sections 2.2 and 2.3 we review the OCMIP-2 formulation of the phosphorus-cycle model. In section 3 we describe our model's solution using the OCMIP-2 reference parameter values. In section 4 we describe the optimization strategy and optimization results. In section 5 we introduce the equilibrium solution sensitivity patterns (S patterns), show how they are computed, and discuss the sensitivity of the model solution to changes in parameter values. In section 6 we examine the sensitivity of the globally integrated new and export production in the form of particulate organic carbon (POC) and dissolved organic carbon (DOC). Finally in section 7 we summarize our results and present conclusions along with suggestions for future developments and applications of our method.

Method of Solution and Model Formulation
[9] After discretization, the governing equations for the biogeochemical fields (section 2.2) can be expressed as a system of differential equations where m T [s, k, a] is the vector of parameters and X is the model state vector consisting of, for example, the concentrations of phosphate and dissolved organic phosphorus at each model grid point. By definition, the tracer distribution is in equilibrium when its time rate-of-change is zero. Finding the equilibrium tracer distribution, X eq , thus reduces to setting dX/dt = 0 in equation (1) and solving Because the source-sink terms are nonlinear functions of X eq , equation (2) is a coupled system of nonlinear equations for which an iterative solution method is needed. We use Newton's method [e.g., Kelly, 2003] and iterate until maxjF(X)j < tol where tol is chosen to be as small as possible given the finite machine precision.
[10] Each Newton iteration requires the inversion of the model's Jacobian matrix, (@F/@X). The Jacobian matrix has a rank on the order of 10 5 but fortunately the matrix is sparse and can be factored efficiently into upper and lower triangular form using a modern sparse solver [Amestoy et al., 2001]. The resulting code is very fast. For example if we initialize the iteration with uniform [PO 4 ] and [DOP] fields, the solver reaches an equilibrium solution with zero time-tendency drift within machine precision in about six to seven iterations, each taking about 5 min on a single processor workstation.

Ocean Transport Model
[11] We use the off-line ocean transport model developed by Primeau [2005]. The model is based on the time-averaged velocity and eddy-diffusivity tensor fields from a dynamical simulation with a full OGCM. The OGCM used is a version of the ocean component of the climate model of the Canadian Centre for Climate Modeling and Analysis, itself based on version 1.3 of the NCAR CSM Ocean Model [Pacanowski et al., 1993]. It uses the KPP vertical mixing scheme [Large et al., 1994] and the GM isopycnal eddy-mixing scheme [Gent and McWilliams, 1990]. The model is forced by monthly momentum, heat and freshwater fluxes obtained from the atmospheric component of the climate model together with restoring of surface temperature and salinity toward their observed climatological values. The off-line transport model (like its parent OGCM) uses second-order centered differences on a grid with 29 levels in the vertical ranging in thickness from 50 m near the surface to 300 m in the deepest level and a $3.75°Â 3.75°horizontal resolution. We tested the numerics of the offline transport model against those of the full OGCM's tracer advection scheme to confirm that they produce the same transport given the same circulation and eddy diffusion tensor. All the biogeochemical simulations presented in the paper use a steady time-averaged transport field with no seasonal cycle.
[12] An important feature of the off-line transport model is that its advection-diffusion operator is available in matrix form. This allows us to compute equilibrium tracer distributions for prescribed source and sink patterns by direct matrix inversion, thus avoiding the need to perform lengthy multithousand year integrations needed for the tracer fields to come into equilibrium with their source and sink patterns. This unique feature of the model has been exploited in previous studies [Primeau, 2005;Primeau and Holzer, 2006] to perform extensive studies of the ventilation and transport properties of the model.

Ocean Biogeochemistry Model
[13] The biogeochemistry model is based on the OCMIP-2 formulation (Najjar and Orr online document) [see also Sarmiento et al., 1988]. The full model has five prognostic variables, phosphate, [PO 4 ], dissolved organic phosphorus, semilabile [DOP], total alkalinity, A T , dissolved inorganic carbon, [DIC], and oxygen, [O 2 ], but in this study we will focus on the phosphorus cycle only.
[14] The governing equations for the phosphorus cycle are where the source-sink terms, S PO 4 and S DOP , are functions of both [PO 4 ] and [DOP] and thus couple the two equations.
[15] The source-sink terms for PO 4 and DOP take the following forms where J PO 4 is organic phosphorus production, z c is the compensation depth, fixed at 75 m, k is the inverse e-folding timescale for the remineralization of dissolved organic matter (DOM), s is the fraction of the production that is recycled as DOM, (1 À s) is the fraction that is exported as particulate organic matter (POM), and a is a constant scaling factor that controls the depth-profile for the remineralization of POM, assumed to be instantaneous.
[16] In order to be able to apply Newton's method to obtain solutions of the model, special care is needed to ensure that the governing equations are differentiable with respect to the dependent variables so that the model's Jacobian matrix is well defined. We thus modified the formulation for the biological uptake of phosphate in new production as follows: where [PO 4 ] obs is the observed phosphate concentration and where l is a small number with the same unit as [PO 4 ]. Note that in the limit of l ! 0, equation (8) reduces to the same formulation proposed by Sarmiento et al. [1988] and implemented in the OCMIP-2 protocol. In practice, we choose l = 10 À5 mol/m 3 , which is large enough for the partial derivatives of the model equations with respect to [PO 4 ] to be numerically well defined, but small enough that there is little difference in the solution obtained with formulations (8) and (9).

Constraining the Total Phosphorus
[17] The ocean transport model is written in flux-divergence form so that it is guaranteed to conserve tracer mass in a timedependent simulation. This conservation property imparts on the transport operator and the model Jacobian a null-space that makes it singular. As a result, the total tracer mass is free to drift from iteration to iteration making it impossible for the Newton solver to converge to a solution with a realistic total mass of phosphorus. An additional constraint equation is thus needed to make the Jacobian matrix nonsingular and invertible. (In practice, round off error makes it possible for our direct solver to invert the ''singular'' matrix, but the resulting total mass of phosphorus is then determined arbitrarily by the round off error.) In a time-dependent model the indeterminacy is removed by prescribing the initial condition which sets the total mass of phosphorus in the ocean. The partitioning of the phosphorus into PO 4 and DOP is then determined as part of the solution. With our iterative solver we can prescribe either the total phosphorus, the total PO 4 , or the total DOP, and allow the solution to determine the total mass of the remaining two tracers. Which constraint is preferable depends on which of the three total masses is best known from observations. We chose to prescribe the total phosphorus as it is easier to interpret the sensitivity of the solution to changes in parameters when total phosphorus is held constant as it would be in a time dependent simulation, but we also implemented the total PO 4 constraint for the parameter optimization study since the total mass of DOP is not well constrained from observations [Karl and Björkman, 2002] and our goal is to constrain the model parameters using the observed [PO 4 ] field. Another option would have been to add the total mass of phosphorus to our list of model parameters to be optimized, but we did not implement this option because of the increase in computational cost associated with performing the optimization in a higher dimensional space.

Equilibrium [PO 4 ] for Control Case With OCMIP-2 Reference Parameters
[18] The equilibrium phosphate distribution for the simulation based on the OCMIP-2 parameter values (Table 1)  [19] The depth averaged [PO 4 ] distribution is contoured in Figure 1. It shows broad agreement between the simulation and the observations: Both have elevated [PO 4 ] levels in the North Pacific and northern Indian Ocean and depleted [PO 4 ] levels in the North Atlantic. However, the total column [PO 4 ] maximum in the equatorial eastern Pacific and in the Bay of Bengal in the Indian Ocean are much higher in the simulation than in the observed field. Furthermore, the observed field shows elevated levels of total column [PO 4 ] in the Gulf of Alaska, the Bering Sea and the Arabian Sea that are absent in the simulation.
[20] Contour plots of the zonally averaged [PO 4 ] for the Atlantic, Pacific and Indian oceans are presented in Figure 2. There is broad agreement between the simulated and observed [PO 4 ] fields. The low- [PO 4 ] waters in the upper thermocline of the subtropical gyres and the bowing upward of the nutricline in the tropics are well captured by the simulation. The tongue of North Atlantic Deep Water (NADW) with low preformed- [PO 4 ] is visible in both the  Figure 2). However the NADW tongue in the model is much thinner than its counterpart in the observations. The observed NADW core is centered at 3000 m whereas the model NADW core is confined above 3000 m by [PO 4 ]-rich Antarctic Bottom Water (AABW) that penetrates too far north in the Atlantic [Primeau, 2005]. This feature is characteristic of coarse resolution models that do not produce upper and lower NADW [Dutay et al., 2002].
Other prominent model deficiencies can be seen in the bottom waters of all the basins where the model overpredicts the [PO 4 ] concentration and in the upper thermocline waters northward of 40°N in the Pacific where the observed [PO 4 ]-nutricline is much sharper and closer to the surface than in the model. The cause of the later is likely due to a misrepresentation of the circulation resulting in upwelling that is too weak in this part of the basin. The results of our optimization study (section 4) will rule out the possibility that these errors are due to the misspecification of the biogeochemical model parameters. The model does however capture the large-scale inter-basin phosphate gradient with high [PO 4 ] in the Pacific and Indian oceans and low [PO 4 ] in the Atlantic.
[21] Overall the simulation captures the broad features of the global phosphate distribution but with some clear discrepancies. By optimizing the biogeochemistry model parameters in the next section we will be able to quantify how much of the unexplained variance can be attributed to misspecification of the biogeochemical model parameters and how much is likely attributable to discrepancies in the circulation.

Parameter Optimization Study
[22] In this section we use global World Atlas 2001 gridded [PO 4 ] [Conkright et al., 2002] to optimize the biogeochemistry model parameters. We do so by adjusting a, s and k so as to minimize the volume integrated squared difference between the simulated and observed [PO 4 ] field, where h[PO 4 ] obs i is the average phosphate concentration of the ocean. The denominator is chosen so that the cost function can be interpreted as the fraction of the spatial PO 4 variance not captured by the model.
[23] Because the effects on the PO 4 distribution of the processes parameterized by the different parameters are not independent, it is difficult to optimize individual parameters separately. In particular, the sensitivity patterns of s and k are highly correlated as we will show in Section 5. We thus optimize all three parameters simultaneously by minimizing the discretized version of (10) using the Nelder-Mead sim- plex method [e.g., Press et al., 1992] as implemented in Matlab's FMINSEARCH.
[24] We performed two optimizations using different constraints: one in which we prescribed the total mass of PO 4 such that the mean PO 4 concentration of the solution is equal to the observed value (2.17 mmol/kg) and another in which we prescribed the total mass of phosphorus such that the mean concentration of phosphorus is equal to its estimated value (2.19 mmol/kg) (see section 2.3). The different constraints result in a small difference in the total amount of phosphorus. As we will show in section 4.2, s and k are not well constrained independently and are thus sensitive to small changes in the total amount of phosphorus. This translates into uncertain estimates for the global inventory of DOP and for the globally integrated new production (section 6). However, the overall model data misfit after optimization, the estimate of globally integrated export production from POC and the optimal value of a are not very sensitive to which constraint is used.
[25] The optimized parameters for each case along with the original OCMIP-2 reference parameters are given in Table 1. The optimized equilibrium [PO 4 ] fields are discussed in the next section followed by a discussion of parameter uncertainty.

Optimal Parameter Set
[26] Our optimized model captures more than 70% of the spatial variance in the observed [PO 4 ] field, a significant improvement over the control run with the OCMIP-2 parameters that captured only 60% of the variance (i.e., an 13% decrease in the RMS error in the PO 4 field).
[27] Most of the improvements (9% of the spatial variance) can be attributed to the combined change in s and 1/k. The small changes in a (Table 1) contribute only an additional 1% in the fraction of variance captured by the model suggesting that the OCMIP-2 reference value was already close to optimum.
[28] The optimal power law exponent, a, for the POP-flux remineralization profile is insensitive to whether the total mass of phosphorus or the total mass of phosphate is constrained during the optimization. Both approaches produce an optimal value of a = À1.0, which is equal to the value first proposed for POC by Suess [1980] and higher than, but still in accord with, the value of À0.9 suggested by Yamanaka and Tajika [1996] and adopted as part of the OCMIP-2 protocol. Yamanaka and Tajika [1996] only tested three values, a = À0.4, À0.9 and À2.0, against the mean-GEOSECS vertical PO 4 profile. Our optimized value is also in accord with the open-ocean-composite value obtained by Martin et al. [1987] from sediment trap data for particulate nitrate, a = À0.988, but smaller than the value for organic carbon, a = À0.858, and also significantly smaller than the mean value of a = À0.68 ± 0.04 obtained for POC by Primeau [2005] in a reanalysis of JGOFS sediment trap data. Apart from the real possibility that the sediment trap data itself is biased, the discrepancy might be due to depth-dependence in the Redfield C:P ratio [Schneider et al., 2003], circulation errors or to the spatial variance of a. Usbeck [1999] used hydrographic data and an inverse approach to estimate a as a function of position for the global ocean and found values that ranged from À0.5 to À2.0 for the remineralization of POC.
[29] Unlike the case of a, the optimized values of s and 1/k are quite sensitive to small changes in the total amount of phosphorus in the model. This, as we shall see in section 5, can be explained by the fact that the sensitivity patterns associated with s and 1/k are highly correlated. For the case where the total phosphorus is prescribed, the result of the optimization is to considerably increase the remineralization timescale for DOP from 1/k = 0.5 years as prescribed by OCMIP-2 to 1/k = 2.8 years. For the case where the total PO 4 is prescribed, the result of the optimization is again to increase the remineralization timescale for DOP, but by a lesser amount to 1/k = 1.0 years. The result of the optimizations on the fraction of production allocated to DOP is to decrease s from 0.67 to 0.61 for the case where the total phosphorus is prescribed but to increase it to 0.74 for the case where the total PO 4 is prescribed.
[30] Remaining model-data misfits are likely attributable to misrepresentation of the ocean circulation as they can not be reduced by further adjusting the biogeochemistry model parameters. However, deficiencies in the formulation of the biogeochemistry model as opposed to a misspecification of the model parameters can not be ruled out. For example, the semidiagnostic form of the new-production parameterization in the model seems to exacerbate circulation deficiencies in the North Pacific subpolar gyre where underestimated surface PO 4 result in a complete shut down of new production which in turn allows nutrients to be exported out of the North Pacific by the more vigorous surface circulation instead of being exported to depth in the form of POP where they would Figure 3. Plot showing the shape of the cost function in the neighborhood of its minimum expressed as a fraction of the total spatial variance that is not captured by the model for the case of total P prescribed with fixed parameters s = 0.61 and 1/k = 2.8 (dashed curve) and for the case of total PO 4 prescribed with fixed parameters s = 0.74 and 1/k = 1.0 years (solid curve).

GB4009
KWON AND PRIMEAU: OPTIMIZATION AND SENSITIVITY STUDY be corralled back into the North Pacific by subsurface currents.

Parameter Uncertainty
[31] Obtaining formal uncertainty estimates on our model parameters requires that the spatial covariance structure of the residuals be known or estimated. This is difficult to do correctly, so we have not attempted to do so, but a sense of the uncertainties can be obtained from the shape of the cost function. For example, Figure 3 shows the shape of the cost function in the neighborhood of its minimum as a function of a. Large values of the cost function indicate unlikely parameter values. Conditional on the model being correct, the curvature of the cost function in the neighborhood of a = À1.0 indicates that the uncertainty associated with this parameter is small, especially considering the large number of degrees of freedom in the PO 4 data.
[32] Figure 4 shows a contour plot of the cost function in the neighborhood its minimum as a function of s and 1/k. The plot shows the contours of C expressed as a fraction of the total variance that is unexplained by the model for the case where the total P is prescribed with a = À1.0 fixed. The shape of the cost function for the case where the total PO 4 is prescribed (not shown) is very similar. Figure 4 shows a broad curving valley defining a family of equally plausible (s, 1/k) pairs. We are therefore able to determine that if the fraction of production allocated to DOP is s = 0.61, as in the OCMIP-2 protocol, then the mean lifetime of this DOP is likely larger than one year and that the simulated PO 4 field is not very sensitive to the exact lifetime of DOP provided it is longer than one year. Alternatively, if the lifetime of DOP is half a year, as in the OCMIP-2 protocol, then the fraction of production allocated to DOP is likely closer to s = 0.85 than to the s = 0.67 value chosen for the OCMIP-2 protocol.

Equilibrium [PO 4 ] for Optimized Parameter Set
[33] The bottom panels of Figures 1 and 2 show the equilibrium phosphate distribution for the optimized parameters obtained using the total PO 4 constraint. The [PO 4 ] field for the total P constraint is rather similar and is not shown.
[34] The most dramatic improvements are in the eastern tropical Pacific (Figure 1) where the optimized parameters reduce the excessive nutrient trapping according to the mechanism first suggested by Najjar et al. [1992] and Anderson and Sarmiento [1995], i.e., by allowing more nutrients to be transported away from the upwelling zone in dissolved organic form. Excessive nutrient trapping in the Bay of Bengal in the northern Indian Ocean is also reduced with the optimized parameters. The optimal parameters do not however completely eliminate the problem confirming the findings of Matear and Holloway [1995] that a proper representation of the circulation is critical to avoiding excess nutrient levels in the eastern tropical Pacific. Aumont et al. [1999] found that increased model resolution improved the nutrient trapping problem but did not completely eliminate it. Together our results suggest that it is the combined misrepresentation of the circulation and misspecification of the biogeochemical parameters controlling the cycling of DOM that are responsible for the excessive nutrient trapping problem.
[35] In the Atlantic (Figure 2), the optimized parameters do little to improve the model-data misfit suggesting that the major discrepancy in the Atlantic can be attributed to errors in the model circulation; the model produces too much Antarctic Bottom Water [Primeau, 2006] and poorly represents the mixing of CFCs and bomb radiocarbon in the far south [Krakauer et al., 2006], a problem that is common in OGCMs of similar resolution [Guilderson et al., 2000;Dutay et al., 2002;Matsumoto et al., 2004]. In the Indian Ocean sector of Figure 2 we see that the optimization produces some improvements, by increasing [PO 4 ] in the southernmost part of the basin and as already mentioned by decreasing the excessive nutrient trapping effect as can be seen in the subsurface [PO 4 ] maximum in the northern part of the basin. However, the location of the maximum is still too shallow in comparison with the observations.
[36] We can partition the fraction of variance captured by the model into high and low latitudes to obtain a quantitative summary of the model improvements. In the latitude band spanning 30°S to 30°N the fraction of variance captured by the model increased from 59% with the reference OCMIP-2 parameters to 72% with the optimized parameters. In the region between 90°S and 30°S it increased from 67% to 75% and in the region between 30°N and 90°N it increased from 58% to 62%. We thus find quantitative improvements in both high and low latitudes, indicating that at least on the largest scales, the optimization does not degrade the fit in one region in order to compensate for circulation errors in another region.

Sensitivity of [PO 4 ] to Changes in Parameters: S Patterns
[37] The difference in tracer concentration at equilibrium for two different values of a parameter, normalized by the Figure 4. Contour plot showing the shape of the cost function as a function of s and 1/k in the neighborhood of its minimum at a = À1.0. The plot shows a fraction of the total spatial variance that is unexplained by the model for the case where the total P is prescribed. The circle indicates the optimal parameter values for the optimization in which total P is prescribed, the square indicates the optimal parameter values for the optimization in which total PO 4 is prescribed.

GB4009
KWON AND PRIMEAU: OPTIMIZATION AND SENSITIVITY STUDY difference in parameter value is a measure of the sensitivity of the equilibrium solution to changes in the parameter. In the limit of an infinitesimal parameter change the resulting pattern is given simply by the derivative of the tracer field with respect to the parameter. We call the resulting sensitivity patterns S patterns and denote them by S a , S s and S 1/k . For example the S a pattern for [PO 4 ] is given by where the subscript eq denotes the fully spun-up equilibrium tracer distribution. S patterns describe the spatial rearrangement of tracers that result from changes in model parameters.
They are particularly useful for elucidating our optimization results.
[38] Our direct solver approach for computing equilibrium tracer fields makes it possible to compute S patterns with little additional computational costs. Through the steady state equation, (2), X eq is defined implicitly as a function of the model parameters. Differentiating (2) with respect to the parameters and solving we obtain where r m (@/@a, @/@s, @/@k). The right hand side, r m F is easily computed analytically from the governing equations. Recall that the model Jacobian matrix, (@F/@X), is available in factored form as part of the Newton solver so that each S pattern can be computed with one additional back-solve -an insignificant increase in computational cost in comparison to the factorization step. (A back-solve is several orders of magnitude faster than factorizing the matrix.) Thus an important advantage to our solution method apart from avoiding lengthy spin-up times is that the sensitivity of the model to its parameters is easily available.
[39] The S patterns for the OCMIP-2 reference parameter case are presented in Figures 5 to 7. The S patterns were computed from equation (12) for the case where the total phosphorus is prescribed and held constant.
[40] In the first row of Figure 5, the zonal-average of S a for [PO 4 ] is shown for each basin. The pattern shows that increasing a results in a slower decrease of the POC flux with depth so that a larger fraction of the flux reaches the deep ocean. The zero sensitivity line, above which [PO 4 ] decreases with increasing a and below which [PO 4 ] increases, varies with latitude and basin but is generally deeper in the Atlantic at $2600 m compared to $1100 m in the other basins. The difference in the Atlantic can be attributed to the formation of North Atlantic Deep Water (NADW) which carries with it water with a reduced amount of preformed [PO 4 ] when a is increased. In the upper ocean, the largest rate of decrease of [PO 4 ] is in the high-productivity low-latitude upwelling regions. The largest rate of increase of [PO 4 ] is in the deep North Pacific consistent with the conveyor-belt picture of the global circulation. The deep water masses of the North Pacific being the oldest show an increased sensitivity to a because they integrate the deep POP fluxes for a longer period of time. However, the highest sensitivities are near the bottom of the North Pacific and Indian oceans whereas the oldest waters for our model are at intermediate depths $1000 -1600 m in the Pacific and $500 -1200 m in the Indian ocean [Primeau, 2005], indicating that not all the sensitivity can be attributed to the accumulation of PO 4 as the water ages. A significant fraction of the sensitivity must be due to the local effect of the deeper remineralization of POP raining down from directly above.
[41] The coupling of biologically mediated transport of PO 4 with the transport by the global circulation results in a net transport of PO 4 from the Atlantic to the Pacific basin. This effect can be seen in the top panel of Figure 6 which shows the column averaged S a for [PO 4 ]. In the Atlantic the net northward flow of PO 4 -depleted water at the surface together with the PO 4 -enriched return flow at depth results in a net export of PO 4 . In the North Pacific and Indian oceans, the overturning circulation is weak, but increasing a pumps more PO 4 into the deep water where currents and eddies (parameterized in the transport model) are weaker than their surface counterparts thus decreasing the efficiency with which the circulation can erode the large-scale PO 4 gradient. The net effect is an accumulation of PO 4 in those basins (Figure 6, top).
[42] The middle and bottom panels of Figures 5 and 6 show S s and S 1/k for [PO 4 ]. Both patterns have very similar zonally averaged and vertically averaged structure. As already mentioned in section 4.1, the similarity explains why the optimization fails to constrain both parameters independently. We show also S a and S 1/k for DOP in Figure 7. The corresponding S pattern for s is very similar to S 1/k and is not shown.
[43] The magnitude of the S patterns for two different parameters cannot be compared directly (they have different ranges and in the case of S 1/k different units also), but a meaningful comparison can be made by multiplying the S patterns by finite parameter differentials. For example, using ds = 0.74 -0.67 and d(1/k) = 1.0-0.5 years, the differences between OCMIP-2 reference parameters and our optimized parameters, we obtain the finite differences  S s ds and S 1/k d(1/k). Each quantity represents a linearized estimate of the contribution of each parameter in producing the difference in the [PO 4 ] field in going from the OCMIP-2 to the optimized parameters. For example, referring to Figure 5 we see that in the equatorial Pacific S s ds yields a maximum reduction in the zonal-averaged [PO 4 ] equatorial region of 0.11 mmol/m 3 and S 1/k d(1/k) yields a zonalaveraged decrease of 0.30 mmol/m 3 in the same region. The combined decreases correspond closely to the actual 0.4 mmol/m 3 difference between the OCMIP-2 reference solution and the optimized solution ( Figure 2).
[44] Increasing either parameter tends to decrease [PO 4 ] over broad regions centered on the equatorial upwelling regions in the eastern part of the equatorial ocean. Increasing s allows a larger fraction of the biological production to be carried away in dissolved form by divergent surface currents as opposed to being exported to depth in particulate form where convergent currents trap the remineralized PO 4 and upwell it back to the surface. Increasing 1/k has a similar effect by allowing DOP to be transported farther from the upwelling region before being remineralized into PO 4 thus decreasing the so-called nutrient trapping effect [Anderson and Sarmiento, 1995;Matear and Holloway, 1995;Aumont et al., 1999]. This effect is also reflected in the S a and S 1/k patterns for DOP shown in Figure 7. Both S patterns are closely related with the pattern of surface biological production, but the S a pattern captures regions of high biological production in terms of regions of negative sensitivity quite faithfully. Increasing a causes POP to be exported deeper in the water column resulting in a local decrease in primary production and local DOP inventory. In contrast the S 1/k pattern which is strongly affected by the wind-driven current systems of the upper ocean shows how increasing the mean lifetime of DOP causes nutrients to be carried farther away from the production regions by lateral currents in dissolved form.
[45] It is interesting to note that S s and S 1/k , the sensitivity patterns associated with parameters that control DOP cycling, are coupled much more closely to the wind-driven circulation than S a , which is associated with POP cycling. The later appears to couple more closely with the global overturning inter-basin circulation. For example, S s and S 1/k have pronounced meridional structure associated with winddriven upwelling and subduction regions ( Figure 5). In contrast S a has a more pronounced vertical contrast associated with the global meridional overturning conveyor-like circulation. This suggests that the S patterns should be particulary useful for interpreting changes in biogeochemical parameters that result form changes in the circulation.

Globally Integrated New Production and Export Production
[46] The globally integrated new production and the globally integrated export production from POC and DOC are important properties of ocean biogeochemical cycles because they allow carbon to be pumped out of the atmosphere and sequestered into the ocean interior. In this section we examine how well each quantity is constrained by the PO 4 observation and our model. Table 2 summarizes the globally integrated new production (NP), export production in the form of POC (EPP) and export production in the form of DOC (EPD), for the OCMIP-2 reference parameters and for the two optimized parameter sets.
[47] As in the work of Gnanadesikan et al. [2002], we define new production as the net uptake of phosphate within euphotic zone. The new production in carbon units is obtained by multiplying the net uptake of phosphate (equation (7)) with the constant Redfield ratio of 117:1 [Anderson and Sarmiento, 1994]. Gnanadesikan et al. [2002] showed that new production is very sensitive to vertical diffusivity in the ocean transport model. Here we find that new production is also very sensitive to the processes associated with the cycling of DOM and POM. As expected, decreasing a causes a shallower remineralization of POM which causes nutrient concentration to increase in the upper thermocline. This results in increased entrainment of nutrients into the euphotic zone thus fueling higher new production. New production is also strongly sensitive to changes in s and 1/k (Figure 8a). Particularly so as s approaches unity and 1/k approaches zero. In this near singular limit a large fraction of new production is exported out of the euphotic zone by eddy-diffusive fluxes ( Figure 8d) and it remineralizes quickly immediately below the euphotic zone where it becomes re-entrained into the euphotic zone to fuel a new cycle of new-production ( Figure 8a).
[48] We find also that our model based estimate of new production constrained by PO 4 observations is very uncertain. The reason being that new production is sensitive to combined changes in s and 1/k for which PO 4 is relatively insensitive. This can be seen most clearly by overlaying a contour of the cost function on top of the plot of globally integrated new production (Figure 8a). Recall that equally plausible parameter values lie along contours of the cost function. Because the contours of the cost function intersect a wide range of new production contours, a wide range of new production values are equally plausible. We can thus conclude that the observed PO 4 distribution poorly constrains new production in our model.
[49] In contrast, the model based estimate of export production (Figure 8b) is well constrained by the optimization procedure because in this case the contours of the cost function are nearly parallel to isolines of export production. Thus the observed PO 4 does provide a strong constraint on the globally integrated export production. Put another way, the global PO 4 is not sensitive to changes in the globally integrated production if these changes do not translate into changes in the export production. For either of our optimized parameter sets, the globally integrated export production is approximately 15 Pg/yr (Table 2).
[50] Figure 8c shows the mean concentration of DOC as a function of the model parameters. The relatively large differences in the predicted mean concentration of DOC for equally plausible model parameters suggest that an accurate estimate of the mean DOP concentration in the ocean would help constrain s and 1/k independently. This would help better constrain the model's estimate of new production. Unfortunately, estimates of the mean amount of semilabile DOP are very uncertain [Karl and Björkman, 2002]. The OCMIP-2 suggested initial DOP concentration is 0.02 mmol/kg. This value would constrain s and 1/k in the bottom of the wedge of plausible parameter values. This is also where the estimates of new production would be most uncertain.

Summary and Conclusions
[51] Our study was motivated by the need to better understand the sensitivity of ocean biogeochemical models to their parameters. To this end we implemented an efficient implicit solver with which to obtain equilibrium biogeochemical tracer distributions without having to perform lengthy time-dependent spin-up integrations. Our implicit solver based on Newton's method is several orders of magnitude faster than the traditional time stepping approach and yields the sensitivity of the equilibrium solution to changes in the model parameters with little additional computational costs. We applied our new methodology to an OCMIP-2 type biogeochemistry model and presented a parameter optimization and sensitivity analysis for the model's phosphorus cycle.
[52] Our optimization study made use of gridded global PO 4 observations to constrain the model parameters. We found that the global PO 4 data provides a useful constraint on a, the exponent in the depth profile of POP flux. The PO 4 data did not however provide a strong independent constraint on s and 1/k, the parameters controlling DOP cycling. Instead, the optimization constrained s and 1/k to a curved one-dimensional space that, interestingly, was roughly parallel to the isolines of export production and perpendicular to the isolines of new production. As a result the model's Globally integrated export production of DOC in PgC/yr. The circle and the square indicate respectively the optimal parameter values for the cases where total P and total PO 4 were prescribed. The wedge-shaped region is the 0.32 contour of the cost function (see Figure 4). Parameters values along isolines of the cost function fit the data equally well. predicted export production is rather well constrained while at the same time its predicted new production is only poorly constrained. Our best estimate of export production as POC was 15 Pg of carbon per year based on an assumed C:P Redfield ratio of 117:1.
[53] The lack of independent constraint for s and 1/k can be understood in terms of the parameter sensitivity patterns (S patterns) we computed for our model. These revealed that large-scale PO 4 sensitivity patterns are very similar for both s and 1/k. The S patterns are also useful for understanding the coupling between biogeochemical fluxes and circulation.
[54] Our optimal values for a, s and 1/k indicate that slightly shallower remineralization depth scale for POP and longer residence time of semilabile DOP at the surface ocean compared to OCMIP-2 parameters improve the explained variance of the observed PO 4 by 10% (9% attributed to the combined change in s and 1/k, and 1% to the change in a). The most dramatic improvements are shown in the eastern tropical Pacific and Indian oceans where the control simulation based on OCMIP-2 parameters greatly overestimated column averaged [PO 4 ].
[55] The remaining deficiencies in the simulated PO 4 fields are likely related to deficiencies in the model circulation. The fact that there is a large remaining fraction of unexplained variance for the model's optimal parameter values highlights the importance of properly simulating physical transport for ocean biogeochemical modeling. In this study we have held the circulation fixed so as to keep the number of parameters we estimate to a minimum. This has the advantage of allowing us to explore how well parameters are constrained by the data. With a different circulation our particular parameter values will no longer be optimal and a new optimization would be needed to adjust them to their new optimal values. Future work should focus on improving the model circulation, including seasonal effects and performing an optimization and sensitivity analysis of parameters associated with the cycling of carbon and alkalinity.