A fast Newton–Krylov solver for seasonally varying global ocean biogeochemistry models

We present a computationally-efﬁcient method for obtaining the fully spun-up state of a seasonally- varying global ocean biogeochemistry model. The solver uses a Newton–Krylov method to ﬁnd the ﬁxed points of the map that assigns to an initial state the value of the model state at the end of a one-year run. Apart from the preconditioner, which we describe in the paper, the method relies on a black-box public-domain Newton–Krylov solver that does not require the explicit construction of the model’s Jacobian matrix. Applied to the PO 4 plus dissolved organic phosphorus (DOP) cycle of an Ocean Carbon Model Intercomparison Project II (OCMIP-2) type model, the solver is more than two orders of magnitude faster than the traditional time-stepping method for spinning up the model. The efﬁciency of the solver is illustrated by using the seasonally varying globally-gridded PO 4 climatology to objectively optimize the parameters that control the mean lifetime of semi-labile DOP and the fraction of new production allocated to DOP. The optimization study demonstrates that the information in the seasonal variations of PO 4 do not provide a signiﬁcantly stronger constraint than the annually averaged data used in previous optimization studies. 2008 Elsevier Ltd. All rights reserved.


Introduction
Ocean-biogeochemistry models have many uncertain parameters that cannot be directly measured in the field or laboratory. These parameters must be constrained by adjusting them so that the model output matches the available observations as well as possible. It is therefore desirable to develop objective parameter optimization methods that go beyond the trial-and-error approach.
The major impediment to applying automatic parameter optimization methods to global ocean-biogeochemistry models is that global models take several thousand years to reach a new equilibrium each time a parameter is perturbed. This timescale for transients to die out and the model to reach its equilibrium is set by the most slowly decaying eigen-mode of the model's advectiondiffusion equation. For a global model the slowest decaying mode has a timescale of roughly 800 years (Primeau, 2005;Khatiwala et al., 2005) and projects directly on the transport of nutrient rich deep water back to the sunlit surface ocean. The long runs needed to reach equilibrium have made it impractical to perform the large number of runs needed to systematically optimize parameters objectively even for relatively coarse resolution models.
Here, we present a new implicit solver that greatly reduces the computational time needed for obtaining equilibrium solutions of a seasonally varying global ocean-biogeochemistry model. We build on the recent work of Primeau (2006, 2008) who introduced the use of an implicit iterative solver based on Newton's method to compute the steady state of an ocean-biogeochemistry model. The key new development we present here is a solver that can obtain cyclo-stationary solutions of a seasonally varying model.
The inclusion of seasonal effects is important for biogeochemistry because the upper ocean exhibits strong seasonal variations in both circulation and export production. Potentially important rectification effects cannot be captured in a steady model. It is known for example, that the T-S properties of the main thermocline reflect more closely those of the end-of-winter mixed layer rather than the annual-mean state and one expects that Stommel's mixedlayer ''Demon" is also at work biasing the coupling between the surface mixed layer and the main thermocline for nutrients, albeit modified by the seasonality of the sinking particle flux (Stommel, 1979;Williams et al., 1995). The first goal of this paper is to present the formulation and implementation of a new solver and demonstrate its efficiency in the context of the phosphorus cycle as implemented in 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). A key advantage of our new solver is that it does not require the explicit construction of the model's Jacobian matrix and should therefore be easy to implement in any ocean biogeochemistry model. In fact the solver we use is a professionally written public domain Newton-Krylov solver written by Kelley (2003). All that is needed to apply the method is to formulate a preconditioner and to implement a simple interface between the model and the black-box Newton-Krylov solver. We give only a basic idea behind the Newton-Krylov algorithm because the implementation details can be found in the book by Kelley (2003). Our focus is on the formulation of the interface between the model and the solver which we describe in Section 2 and the preconditioner which we described in Section 4 after first reviewing the model formulation in Section 3. The formulation of an effective preconditioner is the key to the method's efficiency.
In addition to describing the solver itself, we also illustrate how the fast solver can be used effectively in a parameter optimization study. The optimization study we conduct is motivated by the previous work of Kwon and Primeau (2006, hereinafter referred to as KP06). In their optimization study of the parameters controlling the cycling of phosphorus using the annual-average global PO 4 climatology, KP06 found that the parameters controlling the mean lifetime of dissolved organic phosphorus (DOP) and the fraction of new production allocated to DOP were not well constrained independently. They speculated that perhaps the information in the seasonal variations of the PO 4 climatology that was left out in the cost function for their steady-state model might help to better constrain the parameters. The second goal of the present study is to test this hypothesis. To our knowledge ours is the first such optimization study of a seasonally varying global model. The formulation of the cost function and the optimization results are presented in Section 6.

Equilibrium solutions expressed as the fixed points of a map
In this section, we describe the interface needed to cast the output of a biogeochemistry model in a form that can be used by the Newton solver. For the steady-state case considered by KP06, the equations governing the biogeochemistry model can be written symbolically as where xðtÞ is the model state at time t. The usual approach to spinning-up the model to equilibrium is to integrate Eq. (1) forward it time until the transients die out so that dx=dt ! 0 and x becomes independent of time. To find the equilibrium state using Newton's method one sets dx=dt ¼ 0 and seeks the solution to the resulting time-independent coupled system of algebraic equations, Fðx eq Þ ¼ 0. A sequence of Newton iterates that yields x eq when it converges is where ½oF=ox xn is the Jacobian matrix of partial derivatives of F with respect to the model state x evaluated at x n . States for which the time-tendency are zero within machine precision are typically obtained after 5 or 6 iterations. The efficiency of the solver used by Primeau (2006, 2008) derives from the fact that the Jacobian matrix for finite-difference models is very sparse. Because the matrix has mostly zero elements it can be inverted efficiently using a modern sparse-matrix factorization algorithm. We now wish to apply a similar method to a seasonally varying model. The difficulty is that the equilibrium state is now timedependent so we cannot simply set dx=dt ¼ 0. The governing equation for the seasonally varying model is of the form with an explicit time-dependence in F due to the seasonality of the circulation, temperature field, solar radiation, etc. We are interested in the case where the explicit time dependence of F is periodic with period T ¼ 1 year, i.e. Fðt þ T; xÞ ¼ Fðt; xÞ. In an offline transport model it is easy to ensure that the system satisfies this condition by simply recycling the pre-saved one-year time series of velocity and eddy-diffusivity. In a biogeochemistry model that is run concurrently with ocean dynamics, one can ensure that the system is periodic by recycling a one-year time series of surface fluxes of momentum, fresh-water and heat, and by resetting the dynamical state of the model (temperature, salinity and momentum) to the same state at the beginning of each year. To find the periodic equilibrium solutions of equation (3) we want to recast the problem in a form that can be solved using Newton's method. To this end, we define the map, M : x 0 ! x, that assigns to the initial model state xðt o Þ x 0 , the value of the solution one period later, i.e. xðt o þ TÞ. Evaluating the map, M, is equivalent to running the model forward in time for a one year period starting from the initial state x 0 . The equilibrium solutions of Eq.
(3) correspond to fixed points of M, i.e. Mðx eq Þ ¼ x eq . We note in passing that the usual approach of spinning up an ocean model by timestepping the model forward is equivalent to successive applications of the map M, in other words the model state at the end of one year is used as the initial condition for starting the next year of integration. Spinning-up the model in this way is therefore equivalent to applying Picard's method (e.g. Pozrikidis, 1998), x nþ1 ¼ Mðx n Þ, to find the fixed point of M. An alternative to Picard's method with a better convergence rate is to apply Newton's method to find the roots of Gðx eq Þ Mðx eq Þ À x eq ¼ 0.
The use of a map as in intermediary in the application of Newton's method has been applied successfully in the chemical engineering community to obtain periodic states of cyclically operated chemical processes, (see for example Van Noorden et al., 2003.) Here, we apply the method on a much grander scale to the global ocean viewed as a chemical reactor operated cyclically by the periodic orbit of the Earth around the Sun. Newton's method applied to Gðx eq Þ ¼ 0 yields Unfortunately, the Jacobian of G unlike that of F is not sparse. Even for a low-order finite difference scheme, a tracer perturbation at one grid point can spread to many grid points during the integration of the model over a full year. Consequently, the direct inversion approach used for the steady-state case by Primeau (2006, 2008) is no longer feasible. Instead of explicitly constructing and inverting the Jacobian matrix of G, we use an iterative solver that needs only Jacobian-vector products. These in turn are computed by the solver itself from G using a forwarddifference directional derivative; all that is needed to apply the solver is a subroutine that evaluates G. Evaluating GðxÞ is as simple as running the model forward for one year starting with x as initial condition and then subtracting this initial state from the model state at the end of the one-year run. The nonlinear solver then uses an iterative linear solver, to approximate the solution to Eq. (4). A Newton-Krylov solver therefore consists of two iterative solvers: an outer Newton solver which forms the iterations of Eq. (4) and an inner Krylov solver that uses an iterative approach to ''invert" the Jacobian. The particular linear iterative solver we use is the GMRES algorithm (Saad and Schultz, 1986) as implemented by Kelley (2003). GMRES approximates the solution of a linear system Ax ¼ b with a sum of the form, where r 0 ¼ b À Ax 0 and where x 0 is the initial iterate. The coefficients a k are determined using linear least squares to minimize jjb À Ax n jj. The details of the implementation are taken care by the black-box solver. The convergence rate of the inner iterative solver depends critically on preconditioning the matrix A by pre-multiplying it by an approximate inverse P % A À1 . The closer P is to the exact inverse of A the fewer iterations will be needed. However, there is a trade-off between the computational costs associated with the construction and application of P and the need to minimize the number of inner iterations. The preconditioner we have implemented is the socalled ''physics-based" preconditioner in that it is inspired by the form of the governing equations. We therefore briefly describe the biogeochemistry model equations before returning to the formulation of the preconditioner.

Model description
The OCMIP-2 type phosphorus cycle model is described in detail in KP06. The circulation model is also the same as the one used by KP06. The governing equations of phosphate, PO 4 , and dissolved organic phosphorus DOP are where S PO 4 , and S DOP are the source-sink terms of PO 4 and DOP due to biological uptake and remineralization and where J PO 4 ðu À KrÞ½PO 4 and J DOP ðu À KrÞ½DOP are the advectivediffusive fluxes of PO 4 and DOP for a fluid velocity u and eddy-diffusion tensor K. We have also added a slow phosphate restoring term with PO 4 ¼ 2:2 mmol=m 3 and with s PO4 ¼ 53; 000 years set to the mean residence time of PO 4 in the ocean. This extra sourcesink term eliminates the null space in the Jacobian of G associated with the conservation of total phosphorus in the original model. Apart from making the Jacobian invertible the extra term ensures that the iterative solver will always find a solution whose total phosphate inventory equals the value estimated from observations. Without such a constraint the total amount of phosphorus would be free to drift during the Newton-Krylov iterations. The geological timescale associated with this term ensures that it has little effect on the spatio-temporal PO 4 distribution. The biological source-sink terms are defined as follows in which Here, j PO4 is the organic phosphorus production, also defined as the new production, following Anderson and Sarmiento (1995). It is determined by restoring the modeled ½PO 4 to the observed monthly averaged values, ½PO 4 obs , with a time scale of s ¼ 30 days. The function HðxÞ 1 2 ½1 þ tanhðx=kÞ is a smoothed step-function that is used to shut off production when the modeled ½PO 4 drops below the observed value. This shutting off of the production when nutrient concentrations drop below a threshold introduces the essential nonlinearity in the model. The parameter k controls the steepness of the smoothed-step. In the limit k ! 0; HðxÞ becomes the step-function. We chose a value of k ¼ 10 À5 mmol m À3 which is small enough to give results that are indistinguishable from solutions obtained with the original step-function used in the OCMIP formulation. The parameter j defines the inverse e-folding time scale for the remineralization of DOP. r is the fraction of the production that is allocated to DOP. The remaining fraction, 1 À r is exported out of the euphotic zone as sinking particles (POP). The exponent a controls the attenuation of the POP flux with depth as particles are remineralized and organic phosphorus is returned into inorganic ½PO 4 . The model is solved using an offline model with a one day timestep. The advection-diffusion transport operator uses monthly averaged velocity and eddy-diffusion tensor fields that were obtained from a dynamical OGCM that was spun up for more than 6000 years. Our particular offline model is integrated forward in time using a second order Crank Nicholson scheme for the advection-diffusion terms and a third order Adams Bashforth scheme for the biogeochemical source-sink terms. The off-line model periodically cycles through the 12 u and K fields. It is important to emphasize, that the implementation and efficiency of the Newton-Krylov solver does not depend in any way on the choice of time-integration scheme used in the offline model. We chose this particular time-integration scheme because it was easy for us to implement.
The dynamical OCGM used to generate the circulation fields is a version of the climate model of the Canadian Center for Climate Modeling and Analysis, based on the NCAR CSM Ocean Model (Pacanowski et al., 1993; NCAR CSM Ocean Model Technical Note (NCAR/TN-423+STR) http://www.ccsm.ucar.edu/models/ocnncom/Doc1_3.html). The transport model has 29 vertical levels with thickness of 50 m near the surface to 300 m in the deepest level and a horizontal resolution about 3.75°Â 3.75°. The water mass ventilation properties of the annually averaged circulation are described in detail in Primeau (2005), Primeau and Holzer (2006), Primeau (2006, 2008).

Preconditioner
Our starting point for formulating the preconditioner is based on the insight gained from previous water mass ventilation studies (Primeau, 2005;Primeau and Holzer, 2006) that the long spin-up runs needed for ocean biogeochemistry models to reach equilibrium are due to the slow ventilation timescale for the deep ocean. Even though the preconditioner we use and describe below is ''physics based", it does not have to strictly satisfy any physical or biogeochemical principles. We are free to construct the preconditioner by neglecting most of the biological sources and sinks. We kept only the slow PO 4 phosphate restoring term and the first-order DOP decay rate. These two terms were kept because they eliminate the transport operator's null space. In general, retaining one sink term for each tracer in the model is sufficient to assure that the null space will be eliminated. Because these source sink terms are linear in the model state they do not introduce any state dependence in the preconditioner. In more complex models that include tracers with only a nonlinear sink term one would have to linearize the retained nonlinear sink terms about an appropriate time-independent reference state. Furthermore, since the seasonal variations of the circulation in the deep ocean are not very pronounced, we treated the circulation as if it was steady and used the annually averaged circulation to construct an approximate inverse to oG=ox.
The resulting preconditioner we used is given by P ¼ Q À1 À I with where I is an n Â n identity matrix with n equal to the total number of grid-boxes in the model and where is the advection-diffusion transport operator constructed from the annually averaged velocity and eddy-diffusivity tensor fields. The advection-diffusion equation when discretized using a finite-difference scheme naturally yields a sparse matrix operator because the instantaneous tracer-flux divergence at any given grid point depends only on the tracer concentration in a small neighborhood of grid-cells close to the point at which the flux divergence is being computed. The sparsity of hAi makes it possible to efficiently compute and store the LU factorization of the Q from which the effect of pre-multiplying by P ¼ Q À1 À I can be computed accurately and efficiently (Golub and Van Loan, 1989). Some other important properties of P are that: (1) it is independent of the biogeochemistry model state as well as time independent and so needs to be constructed only once, (2) the numerical scheme used to construct hAi need not match exactly the numerical scheme used for the time dependent tracer transport model so that the time-dependent scheme could in principle use a high-order upwind scheme with flux limiters, and (3) the matrix Q has a block triangular structure which means that only the diagonal blocks of the form hAi À D, where D is a diagonal matrix whose coefficients are the linear Taylor series coefficients of the retained sink terms, need to be factored in order to apply P. This last point suggest that when constructing Q for biogeochemical models that carry a large number of tracers it should be especially advantageous to neglect source-sink terms so as to retain a block-triangular structure. The above preconditioner is applied by pre-multiplying G by P. This is done by solving QG ¼ G using the precomputed and stored LU decomposition of Q and passingG À G instead of G to the blackbox Newton-Krylov solver. Left-multiplying by P does not change the fixed points of M but greatly accelerates the convergence rate of the solver as we will see in the next section.

Convergence rate to equilibrium
The convergence rate of the new solver is shown in Fig. 1 where the root mean square (rms) drift of the model solution is plotted as a function of the number of evaluations of the map M, i.e. as a function of the number of one-year runs used. The figure compares the convergence rate of the traditional time-stepping method (Picard iterations) which takes in excess of 3000 years to reduce the rms drift below 10 À9 mmol m À3 yr À1 , the Newton-Krylov solver without the application of the preconditioner which requires approxi-mately 500 years of integration to reduce the drift below 10 À9 mmol/yr and the Newton-Krylov solver with the application of the preconditioner described in Section 4 which reduces the rms drift to below 10 À9 after only 24 years of integration. Since the computational cost associated with the application of the preconditioner and the overhead associated with the Newton-Krylov solver is negligible compared to a one-year run of the model the preconditioned Newton-Krylov solver is more than two orders of magnitude faster than the traditional time-stepping method for spinning up the biogeochemistry model. All three cases shown in Fig. 1 were initialized with the same observed climatological PO 4 distribution and a uniform DOP distribution. If however one has a good initial iterate, the average number of years of integration needed for the preconditioned Newton-Krylov solver to converge can be substantially less. Good initial iterates are typically available when one uses the solver with a parameter optimization routine. This will be illustrated in the next section.
Depth versus time plots of the fully equilibrated DOP and PO 4 concentrations for the model with parameter values (r ¼ 0:74, 1=j ¼ 1 yr;a ¼ À1:0) are shown for the N Atlantic and Eastern Equatorial Pacific in Figs. 2 and 3. In the Eastern Equatorial Pacific the solution shows only a weak seasonal cycle and surface waters have elevated DOP levels throughout the year. The new-production for the model solution is very high and does not vary much with the seasons in the Eastern Equatorial Pacific ( Fig. 4a and b). The over estimated production in the upwelling region of the eastern equatorial pacific is a common problem for coarse resolution models and is believed to be associated with errors in the model circulation (Najjar et al., 2007). In the N Atlantic, the solution has a strong non-sinusoidal seasonal cycle with a mixed layer that penetrates down to approximately 400 m. At depth, the solution shows only a weak seasonal cycle due to the strong attenuation of the particle remineralization profile. In the surface the DOP concentration is anti-correlated with the PO 4 concentration as expected form the periodic mixing of deep waters that have high PO 4 and low DOP levels and from the fact that DOP the production of DOP is associated with the consumption of PO 4 . We note that the model solution is fully equilibrated even in the deep ocean and shows no drift associated with the slow adjustment from the initial condition.

Optimization of the parameters controlling the cycling of DOP
In this section, we illustrate the usefulness of the Newton-Krylov solver in an optimization study of the parameters that control the cycling of DOP. The parameters in question are r the fraction of new production allocated to DOP and j the first-order e-folding remineralization rate of DOP. When these parameters were optimized by minimizing a cost function formed from the steady model and the annually averaged PO 4 observations, KP06 found that these two parameters were strongly correlated and could not be constrained independently. Our goal here is to test the hypothesis that the information in the seasonal variations in the PO 4 datainformation that was ignored in the steady-state cost function used by KP06 -can help to better resolve the parameters.
To test this hypothesis we form a cost function that measures the distance of the monthly averaged model state from the monthly averaged PO 4 climatology from the World Ocean Atlas 05, (Garcia et al., 2006). The shape of the cost function near its minimum determines how well the parameters are constrained. If the information in the seasonal variations helps to better constrain the parameters, including seasonal variations in the cost function should either increase the curvature of the cost function near its minimum, or decrease the amount of parameter interactions by log [ rms drift(mmol m-3 year-1)] 10 Fig. 1. Root mean square differences of xðt þ TÞ À xðtÞ, where T ¼ 1 yr, for the usual time-stepping spin up (solid line) and for the Newton-Krylov solver with no a preconditioner (open circles) and for the Newton-Krylov solver with the pre-conditioner (solid squares). Although it is difficult to determine from the figure, a total number of 24 simulation years were needed for the pre-conditioned Newton-Krylov solver to converge. aligning the axes of the elliptical bowl around the minimum with the parameter axes.
The cost we use is where DV i is the volume of the ith model grid box, ½PO obs 4 i;k is the observed PO 4 concentration for the kth month interpolated to the ith model grid point and ½PO 4 i;k ðr; jÞ is the simulated equilibrium PO 4 concentration in the ith grid box averaged over the kth month (the model uses a one day time-step). The cost function is scaled by the spatio-temporal variance of the observed field where h½PO obs 4 i ¼ 1 12 is the spatially and annually averaged PO 4 concentration.
To evaluate the cost function we computed equilibrium solutions on a grid of r-j parameter values using a time-step of dt ¼ 1=72 years. For this computation the equilibrium solution at one set of parameter values was used as the initial iterate for a nearby set of parameter values. The same preconditioner was used to compute all the solutions. Fig.5a shows the number of simulation years needed by the solver to reach equilibrium as a function of r and j. The plot shows that with a good initial iterate the typical number of years the model must be run to reach equilibrium is between 40 and 15 years. If one applied an automatic minimization routine to find the optimal parameter values most of the needed equilibrium solutions would be computed in the neighborhood of the minimum of the cost function where for our example the number of years of simulation needed is around 20, making the Newton-Krylov solver more than 150 times faster than the straight time-stepping approach.
The resulting cost function is plotted as a function of r and j in Fig. 5b. The contours indicate the fraction of the spatio-temporal variance that was not captured by the model. The plot shows that the optimal value for r is between 0.4 and 0.75, and that 1/j is not well constrained. The model captures slightly more than 64% of the spatio-temporal variance. This fraction is slightly less than the fraction (70%) of the spatial variance alone captured by the steady model (KP06). The shape of the cost function is also very similar for the steady and seasonally varying model. The cost function for the steady-state model is shown in Fig. 4 of KP06. The conclusion must therefore be that in the context of the OCMIP-2 type model, the information in the seasonal cycle of the PO 4 data does not provide a stronger constraint than the annually averaged data. This result emphasizes the importance of obtaining actual DOP or DOM data so that j and r can be properly constrained.

Discussion and summary
Global ocean biogeochemistry models typically need long spinup runs of more than 3000 years before the deep ocean reaches equilibrium with the surface. This long spin-up time is due to the slowly decaying modes of the circulation model's advective-diffusive operator. In order to decrease the prohibitive computational costs associated with the need to repeatedly re-equilibrate the model each time a parameter is perturbed in the process of optimizing parameters, we have developed an efficient Newton-Krylov solver that is up to two orders of magnitude faster than straight time-stepping. The solver we have developed uses a professionally written public-domain Newton-Krylov solver and does not require any modification to the existing ocean biogeochemistry model. In order to apply the method all that is needed is an interface that forms the difference between the model state at the beginning and the end of a one year run and then solves a linear system of equations in order to apply a preconditioner. Solving the linear system of equations effectively inverts the circulation model's annually averaged transport operator and thus provides an implicit treatment of the model's slow advective-diffusive transport modes. The applicability of implicit solvers for reducing the computational costs of spinning up global biogeochemistry models initiated by Primeau (2006, 2008) is greatly increased by our new solver because it yields equilibrium solutions to seasonally varying models.
In this study, we have explicitly constructed a preconditioner for the phosphorus cycle of an OCMIP-2 type ocean biogeochemistry model with two active tracers. When we started this project, our search in the literature on Krylov methods yielded little in terms of general guidance or specific examples for constructing ''physics based" preconditioners for advective-diffusive systems. We hope that by reporting a specific example of an effective preconditioner for an ocean transport model we might stimulate further research on general methods for constructing efficient preconditioners. Our results here suggest that the preconditioner does not necessarily need to depend on any of the complex biogeochemical parameterizations in order to be effective. Furthermore, our proposed preconditioner is based on the inverse of a block-triangular matrix. This means that only the diagonal blocks of the form hAi À D, where D is a diagonal matrix, need to be factored. This fact means that the computer memory and cpu time needed to form, factor, store and apply the preconditioner will scale linearly with the number of tracers in the biogeochemical model. We therefore expect that preconditioners such as ours, that are based on the annually averaged advection-diffusion transport operator, will continue to be effective in more complex biogeochemistry models as long as the biogeochemical dynamics of these models does not involve any intrinsic timescales comparable to or longer than the longest timescales associated with the circulation, i.e. $1000 years. Whether or not it becomes necessary to modify the preconditioner should the intrinsic dynamics of the biogeochemistry model include such long timescales is difficult to predict, but we expect that the Newton-Krylov solver will perform as least as well as the un-preconditioned case which achieves a non-trivial factor of 10 speedup compared to straight time-stepping.
We have also tested that the LU decomposition of blocks of the form hAi À D with the same sparsity pattern as would be obtained from our model with a 1 Â 1 resolution with 24 levels can be computed using less than 32 GB of memory in under 1.5 h using one Intel(R) Xeon(R) 2.33 GHz cpu. Furthermore, this factorization needs to be computed only once. The backsolves needed when applying the preconditioner at the end of each year of simulation took less than 10 s. The cost of computing and applying the preconditioner is therefore substantially less than running the model for one year. These results suggest that our solution method is expected to substantially speed up models with a resolution comparable to what is currently used for global climate models.
Our solver should also be effective for finding solutions to models that include periods other than the seasonal cycle. In particular, our solver should be applicable with little or no modification to a model that includes a diurnal cycle in addition to the annual cycle. If leap years are not taken into account, the solver should apply without modification. If leap years are included then one would have to use the solver to find the fixed point of the map that propagates the model state forward for four years instead of just one. The effects of longer period phenomena that repeat after an integral number of years could also be included by applying the solver to the map that propagates the model state forward for n years where n is the number of years in the long period cycle. This would make it possible to include some representation of the effects of low-frequency variability on biogeochemical cycling. One could, for example, include a periodic representation of the four year component of El Niño Southern Oscillation (ENSO).
To illustrate the new solver's efficiency we have used it in a parameter optimization study of the parameters controlling the cycling of DOP. Previous studies aimed at constraining r, the fraction of new production allocated to DOP and 1/j and the mean lifetime of semi-labile DOP, have been conducted in the context of steady models without a seasonal cycle (Yamanaka and Tajika, 1997;Primeau, 2006, 2008). In the present study, we have explicitly computed a cost-function that measures the difference between the model simulated PO 4 and the monthly PO 4 climatology from the World Ocean Atlas (2005). By doing so, we were able to demonstrate that the information in the seasonal variations of the PO 4 observations do not significantly improve the resolution of r and j. This highlights the need for direct DOP and DOM observations in order to constrain these parameters.
In this study, we did not make use of an automatic minimization algorithm to optimize the model parameters, but there are no barriers to doing so. Automatic minimization algorithms can be used synergistically to gain additional computational savings. By using a lower tolerance on the Newton-Krylov solver in the early stages of the minimization and increasing the tolerance only  5. (a) The number of simulation years needed for the pre-conditioned Newton-Krylov solver to reach equilibrium as a function of r and j. The same preconditioner was used to compute all the solutions. (b) Cost function expressed as the fraction of unexplained spatio-temporal variance in the observed PO 4 field as a function of r and 1=j the parameters that control respectively the mean life-time of semi-labile DOP and the fraction of new production allocated to DOP as opposed to sinking particles.
in the final stages the overall number of simulation years needed to find the optimal parameters can be further reduced. In the future we plan on using the new solver in combination with an automatic minimization routine to objectively optimize the large and increasing number of parameters in state-of-the-art ocean biogeochemistry models.