An offline implicit solver for simulating prebomb radiocarbon

It takes several thousand years for the deep-ocean concentration of natural radiocarbon to come to equilibrium with surface ﬂuxes, making it computationally too expensive to routinely simulate it with moderate-to high-resolution ocean models. We present an implicit solver for computing prebomb D 14 C that requires the equivalent of only a few tens of model years to reach equilibrium. The solver uses a Newton–Krylov algorithm with a preconditioner based on a coarse-grained annually-averaged tracer-transport operator. Coarse-graining provides a general approach for developing preconditioners for models of increasing resolution. We implemented and tested the solver for the ocean component of the Community Earth System Model (CESM) with a nominal horizontal resolution of 1 (cid:2) (cid:3) 1 (cid:2) and with 60 vertical levels. Simulated D 14 C values are in good agreement with observations at the surface and in the North Atlantic, but the deep North Paciﬁc simulated values show a substantial bias, with prebomb radiocarbon D 14 C values translating to ages that are twice the observationally based estimate. This bias is substantially larger than published simulations obtained with coarser resolution models, suggesting that increasing model resolution does not automatically improve the ﬁdelity of the deep ocean ventilation processes. We therefore recommend that natural D 14 C be used as a deep-ocean ventilation metric for critically evaluating deep ocean circulation.


Introduction
Assessing the fidelity of the ocean circulation in climate models is an important model development step.Models which are used to study climate trends and the evolution of the ocean's CO 2 uptake in a warming world need to be adequately tested before their results are used to influence science and national policy.Because the ocean interacts with the rest of the Earth system over a wide range of spatial and temporal scales there is no unique metric by which to judge the quality of the simulated circulation, but at long timescales the role of the deep ocean becomes increasingly important because of its capacity to store vast amounts of heat and carbon.To assess the rate at which the deep ocean communicates with the surface ocean and the atmosphere, ocean modelers have long recognized the utility of simulating natural radiocarbon (see for example the references listed in Tables 1-3).The availability of globally-gridded natural radiocarbon observations, (GLODAP, Key et al., 2004), has made radiocarbon simulations especially useful for identifying biases in the ventilation of the deep ocean.
For example, Doney et al. (2004) evaluated and contrasted radiocarbon simulations done using 13 different models as part of the Ocean Carbon Model Intercomparison Project Phase 2 (OCMIP-2) and found that errors in the simulated radiocarbon could be attributed to biases in the circulation, and that a significant part of the differences among the models could be tied to differences in sub-gridscale parameterizations.Duffy et al. (1997) and England and Rahmstorf (1999) used D 14 C simulations to evaluate the effect of the Gent-McWilliams (GM) parameterization on the ventilation of the deep ocean, and found that it tended to limit the depth of convection at high latitudes.Further confirming this result, Gruber et al. (2001) found using D 14 C tracer simulations, that the GM parameterization tended to give a sluggish deep circulation in their model.They showed that a 6-fold increase in the vertical diffusivity south of 50°S was needed to reduce excessive stratification and thereby improve the ventilation of the deep Pacific.Gnanadesikan et al. (2004) used D 14 C simulations to evaluate the sensitivity of variations in the vertical and horizontal diffusivities in the Modular Ocean Model version 3 (MOM3).They found that both lateral and vertical mixing processes can affect the resulting D 14 C distribution, but also that changes in the surface forcing can have a major impact.Butzin et al. (2005) found a strong influence of Antarctic sea ice formation on the circulation, based on D 14 C simulations.
Despite its utility, radiocarbon simulations are not routinely done by climate model developers.Simulating natural radiocarbon is a considerable computational challenge because of the long timescales with which deep ocean radiocarbon equilibrates with the atmosphere -some models that participated in the OCMIP-2 study reported tracer equilibration times for radiocarbon of more than 4000 years (Orr, 2002).Furthermore, the computational challenge increases dramatically with increasing resolution because of Courant-Fredrichs-Lewy (CFL) stability criterion restrictions on the model time-step size.As far as we know, none of the ocean components in the current suite of IPCC-class climate models with a resolution of 1 Â 1 or higher have simulated natural radiocarbon to help calibrate the deep-ocean ventilation rate in their models.
Here we present a fast offline implicit solver for simulating natural radiocarbon in the ocean component of the Community Earth System Model (CESM).Our solver uses a Newton Krylov method and preconditioning strategy that are similar to the ones first suggested by Li and Primeau (2008) for biogeochemical tracers, but adapted for the increased memory requirements associated with the higher resolution of the CESM ocean model component, i.e. the Parallel Ocean Program version 2 (POP2), with a nominal 1 Â 1 horizontal resolution and 60 vertical levels (Smith et al., 2010).The cyclo-stationary Newton-Raphson solver and the preconditioner for the iterative Krylov-subspace linear-system solver are described in Section 6.The offline tracer-transport model and the preconditioner are based on monthly and annually averaged tracer-transport matrices, which are constructed (see Section 4), using output from a POP2 simulation, which we refer to as the parent model (Section 3).Using the new implicit solver we are able to compute the natural prebomb-D 14 C equilibrium distribution by running the offline model through only 23 annual periods (Section 7.1).The formulation of the natural radiocarbon model is presented in Section 2 and the time integration scheme used by the offline model to simulate a seasonal cycle is presented in Section 5. Section 7.2 discusses the impact of neglecting the seasonal variations in ocean circulation on the simulated D 14 C distribution.
We compare our D 14 C simulations to the GLODAP observationally based estimates in Section 7.3, providing insight into the behavior and biases of the parent model.We also compare our results to previous radiocarbon modeling studies in Section 8 before presenting our conclusions in Section 9.

Governing equations for 14 C
Following Toggweiler et al. (1989a) we formulate our prebomb radiocarbon model in terms of the ratio R ¼ 14 C= 12 C expressed in arbitrary units in which the prebomb atmosphere is set to where u is the fluid velocity vector, K is the eddy-diffusivity tensor, and is the radiocarbon source-sink function including both radioactive decay with rate constant k ¼ ð8266:6 yearsÞ À1 and air-sea fluxes  parameterized as a source in the top layer of the model.In the airsea flux parameterization l ¼ ð2 yearsÞ À1 and Dz 1 ¼ 10 m.This simplified model for the air-sea exchange of radiocarbon corresponds to the one used in Experiment A of the study by Toggweiler et al. (1989a).It neglects any spatial inhomogeneities in the air-sea flux of radiocarbon due to variations in the air-sea flux of CO 2 .See Section 7.3 for further discussion on the adequacy of this simplified parameterization.
After having been discretized, the governing equations can be expressed as a system of ordinary differential equations in matrix form as follows: ð3Þ in which M is a time dependent matrix defined as where TðtÞ is the advection-diffusion tracer transport operator expressed in matrix form and l is a diagonal matrix whose elements are equal to l for grid-cells in the surface layer of the model and 0 otherwise.Note that the air-sea exchange term has been split into the term lR, which is incorporated into the M term, and the constant term lR atm .
The arbitrary units of R can be converted to the standard D 14 C notation as follows: and R can also be converted to a radiocarbon age as follows: With a half-life of 5730 years and a relatively constant 14 C= 12 C ratio during the Holocene, natural radiocarbon expressed as D 14 C is a useful tracer for quantifying the ventilation of the deep ocean.

Parent ocean model: CESM POP2
The ''parent'' ocean model from which we constructed the seasonally varying tracer-transport operator is based on a prognostic simulation of the ocean component, POP2, of the Community Earth System Model (CESM) (Smith et al., 2010).The POP2 configuration we used has a dipolar grid with the North Pole displaced into Greenland, with the transition from the Mercator grid starting at the Equator (Smith et al., 2010).The computational mesh has N x ¼ 320 and N y ¼ 384 grid points in the nominally eastward and northward directions and N z ¼ 60 grid points in the vertical direction, for a total of 4,241,988 wet grid points.The vertical resolution ranges from 10 m for the top 10 layers of the model and increases to 250 m near the bottom.The thickness of the top-most layer is allowed to vary to include regional and temporal variations in sea surface height.
The dynamical model was ''spun-up'' from rest for 150 years using the Common Ocean-ice Reference Experients (CORE) climatological forcing (Griffies et al., 2009;Large and Yeager, 2009), in which the same seasonal cycle is repeated every model year.The relatively short dynamical spin-up, of only 150 years, was chosen as a compromise between the availability of computational resources and the time needed for the transients in the momentum equations to decay.Although the model is still drifting after 150 years, the relative change in the magnitude of the major currents has decreased substantially and is typically less than 0:1% per year at the end of the spin-up.After the spin-up, the dynamical ocean model was run for an additional year during which we saved all quantities (see Section 4) needed to construct the model's advection-diffusion transport operator in matrix form with a monthly time resolution.This procedure allows us to capture the seasonal cycle, which accounts for the largest part of the transport variability in the ocean.We are thus assuming that interannual variability has a second order effect on the transport of tracers in the ocean, but the error associated with this assumption remains to be quantified.Although we used only one year of the OGCM circulation to capture the transport operator, it is feasible to use an average over several years; this would average out some of the interannual variability.We hope to explore this possibility in the future.
The transport of radiocarbon and other tracers in the model is achieved by a combination of explicitly resolved currents and by a suite of parameterized sub-gridscale transport processes, which include spatially-varying isopycnal and diapycnal diffusion, eddyinduced advective currents, and non-local transport by unresolved density currents over a few key sills where deep waters are formed.Recently incorporated changes in POP2 (Danabasoglu et al., 2012) include a sill overflow parameterization (Briegleb et al., 2010;Danabasoglu et al., 2010), a near-surface eddy flux parameterization (Danabasoglu et al., 2008), a prescription for lateral tracer diffusivities that vary in the vertical (Danabasoglu and Marshall, 2007), a sub-mesoscale mixed layer eddy parameterization (Fox-Kemper et al., 2011), incorporation of a parameterization for abyssal tidal mixing (Jayne, 2009), and a new method for determining background vertical diffusivities (Jochum, 2009).
The model uses the Gent and McWilliams (GM) (Gent and McWilliams, 1990) isopycnal transport parameterization in its skew-flux form with a diffusion tensor that is rotated to allow for different along-and across-isopycnal diffusivity coefficients.The diffusivity has a vertical dependence that varies with the stratification, which can lead to large values in the upper ocean, but much smaller values in the deep (Danabasoglu and Marshall, 2007).The effects of diabatic mesoscale fluxes within the surface are taken into account based on the near-boundary eddy flux parameterization of Ferrari et al. (2008) as implemented in Danabasoglu et al. (2008).High isopycnal diffusivity values, here as large as 3000 m 2 s À1 , in the upper ocean represent the effects of the vertical divergence of eddy stress (Danabasoglu and Marshall, 2007).The values are diminished to 300 m 2 s À1 by a depth of 2000 m.Isopycnal diffusivity coefficients are tapered for isopycnal slopes greater than 0.3 in the interior ocean.The thickness and isopycnal diffusivity coefficients vary together vertically (Danabasoglu et al., 2012).
Vertical diffusivities are determined dynamically using the Kprofile parameterization (KPP) (Jochum, 2009;Danabasoglu et al., 2006;Large et al., 1994).With KPP, the transport of tracers by unresolved convective events is represented by increasing the vertical diffusivity coefficient to 1 m 2 s À1 where the water column becomes statically unstable.The KPP parameterization also includes a background vertical diffusivity due to internal waves that varies with latitude.It is symmetric about the Equator with a minimum value of 0:01 Â 10 À4 m 2 s À1 at the Equator, rises sharply to its mode value of 0:17 Â 10 À4 m 2 s À1 , but with a sharp peak to its maximum value of 0:30 Â 10 À4 m 2 s À1 near 30 of latitude.This background vertical diffusivity does not vary with depth, but the parameterization does include enhanced vertical diffusivities over rough topography with a maximum value of 100 Â 10 À4 m 2 s À1 to mimic the effect of abyssal tidal mixing.
Fig. 1 shows the resulting vertical diffusivity parameter, K V , for the December through March (DJFM) season and the June through September (JJAS) season.The seasonal changes in K V are strongest in the North Atlantic where there is winter-time deep convection.The Southern Ocean shows much less change with the seasons.
Fig. 1 also shows the zonally averaged isopycnal diffusivity coefficient, K I , for the DJFM and JJAS seasons.Generally, the values are lowest along the bottom near land masses and are largest in the seasonal thermocline.There is a substantial amount of seasonality in the value of the K I coefficient.In the North Atlantic, the value of K I rises to over 1000 m 2 s À1 in winter with enhanced values extending down to $4000 m.Southern Ocean K I values are largest during austral winter, but are generally lower than their counterpart in the north.
Also included in the model is a sill overflow parameterization that attempts to capture the unresolved transport of tracers by narrow density currents along the bottom over four key sills: in the Denmark Strait and the Faroe Bank Channel in the North Atlantic, and in the Ross Sea and Weddell Sea in the Southern Ocean (Briegleb et al., 2010;Danabasoglu et al., 2010).The parameterization produces ''non-local'' transport from a ''source'' region upstream of the sill and from an ''entrainment'' region partway down the slope to a ''product'' region in the abyssal basins downstream of the sill.The amount of water flowing from the source to the entrainment and product regions, and from the entrainment to the product region, is determined dynamically depending on the density of waters in the source, entrainment, and product regions.For our particular simulation the surface forcing resulted in watercolumn density profiles that produced sill overflow only in the North Atlantic sills.The conditions in the model were such that there was no sill overflow in the Southern Ocean.
To give a sense of the advective transport responsible for ventilating the deep ocean, Fig. 2 shows the annual, DJFM, and JJAS averages for the meridional overturning circulation (MOC), for the global ocean and separately for the Atlantic basin.The Atlantic basin characterized here includes the entire North Atlantic region: the Labrador Sea, the Greenland-Iceland-Norwegian (GIN) Sea, and the Arctic in addition to the Atlantic itself.The Atlantic MOC shows a northern sinking branch that penetrates down to approximately 2750 m and carries up to $20 Sv.The global MOC shows a strong seasonal variability at low-latitudes with a peak-to-peak amplitude of $100 Sv that is not evident in the Atlantic MOC, indicating that the this variability is confined to the Indo-Pacific basin.This variability, likely due to the changes in the surface wind associated with the Indian Monsoon, leads to a strong tropical MOC cell that changes directions with the seasons and vanishes in the annual average.

Tracer-transport matrix
To capture the seasonal variability of the explicitly resolved advective currents and parameterized eddy-diffusive processes we constructed monthly-averaged tracer-transport matrices, which were formed by separating the transport matrix into three parts, where A includes the effect of advection and the tracer-transport effects of the sill overflow parameterization, D includes the effect of vertical diffusion, and H includes the effect of horizontal diffusion.The sign convention is the same as in the POP manual (Smith et al., 2010).The separation was done primarily to facilitate the reconstruction of the operator using impulse response functions, but it also makes it possible to run the offline model with the same computationally efficient time-stepping scheme as in the parent model, which is based on an implicit scheme for D and different explicit time-stepping schemes for H and A (see Section 5).The D matrix was hand coded using the same second-order centered-difference scheme used in the parent model, as was done in developing the offline tracer-transport matrix of Primeau (2005).The vertical diffusivity coefficients were obtained by time-averaging the instantaneous vertical diffusivity values determined dynamically during a one-year run of the parent model.
During the same run of the parent model we also computed a set of impulse response functions (IRFs) which we used to reconstruct A and H following an approach similar to that described in Khatiwala et al. (2004).Because our parent model used a linear advection scheme without any flux limiters we did not need to use vertically smoothed impulses as recommended in Khatiwala et al. (2004).Furthermore, because the advection scheme in the version of POP2 that we used is based on a third-order upwind spatial discretization with an explicit time-stepping scheme, the response due to A of a passive tracer impulse in grid box ði 0 ; j 0 ; k 0 Þ is guaranteed to extend no further than to two adjacent grid boxes in each direction after a single time-step.(In contrast, the action of D, which is time-stepped using an implicit scheme in the parent model, spreads the response to all 60 levels in the vertical in a single time-step.)The limited spread of the IRF makes it possible to compute the response to impulses separated by four grid cells in any direction simultaneously using a single tracer without any interference among the responses.the time-tendency due to an impulse initial condition which we reinitialize at the beginning of each time step to Cði; j; kji 0 ; j 0 ; k 0 Þ ¼ Maskði; j; kÞd i 0 ;iðmod 5Þ d j 0 ;jðmod 5Þ d k 0 ;kðmod 5Þ ; in which d ij is the Kronecker delta function and Maskði; j; kÞ is a mask equal to one everywhere except for the ''source'' and ''entrainment'' regions of the sill-overflow parameterization.Mask is used to zeroout C for impulses in the source and entrainment regions where the sill-overflow parameterization has been implemented.This was necessary to eliminate the interference in the response produced by the non-local transport into the product areas.(Two separate set of tracers were used to construct the IRFs for each grid box in the source and entrainment regions of the four sills.These additional IRFs were used to reconstruct the sill-overflow effect of A.) In practice, memory constraints and the particular implementation of the IRF module required us to repeat the one-year run seven times to collect all the information needed to construct the twelve monthly transport matrices.The resulting average tendencies, C 0 ði; j; kji o ; j 0 ; k o Þ, for each of the 125 tracer impulses could then be separated into the individual IRFs, using the mask function, Iði; j; kji 0 ; j 0 ; k 0 Þ 1 if ðji À i 0 j 6 2 or ji À i 0 j P N x À 2Þ; jj À j 0 j 6 2; and jk À k 0 j 6 2; The IRFs were then reorganized into the elements of A with the column index determined from ði 0 ; j 0 ; k 0 Þ and the row index determined from ði; j; kÞ.A similar process with a separate set of IRFs was used to reconstruct H.

Enforcing local mass-conservation
By assuming that the year-to-year changes in circulation are negligible to first order we can use the same tracer-transport matrices computed for one particular year to efficiently compute the spin-up of passive tracers without having to solve the momentum and continuity equations.One difficulty with this assumption is that the sea-surface height displacements do not average to zero in any given year.In the dynamical model strong dynamical feedbacks prevent the development of persistent trends in the sea-surface height, but in our offline transport model with a fixed circulation captured from one particular year the local divergence of the flow field that causes the sea-surface height to move up or down does not cancel if we recycle the same circulation year after year.This divergence, although small, accumulates to produce unphysical results for an implicit solver, which produces the solution that would be obtained if the tracer model was run for an infinitely long time.
To eliminate this problem we have modified the horizontal advective flow field in the surface layer of the offline model to ensure that the annual-mean vertical velocity vanishes at the sea surface.(Only gridboxes in the surface layer change volume due to variations in the sea surface height.Consequently velocity corrections are required only for the horizontal velocity in the top layer of the model.)We obtained the correction to the velocity field in terms of a velocity potential, ðũ; ṽÞ ¼ rU, with the potential determined by solving the two-dimensional Poisson equation where gðtÞ is the sea surface displacement and the overline indicates an annual average.Eq. ( 11) was solved subject to no-flow boundary conditions in the domain defined by the basin geometry of the top layer of the model.The resulting velocity corrections were generally small compared with the monthly surface current speeds -the rms change in current speed is $3%.

Annual average circulation
The annually-averaged transport operator, T, is computed by averaging the transport operators constructed from the monthlyaveraged circulation fields for each of the components, A; H, and D, which are added together, as shown in Eq. ( 7) to give the annually-averaged steady-state transport operator, T. The resulting annually-averaged transport operator is then used in Eqs. ( 4) and (3).The steady-state solution of the resulting equation is obtained by setting the time-derivative to zero and directly inverting the time-independent M using the MUMPS multifrontal sparse-matrix solver (Amestoy et al., 2001).The 14 C results obtained using the annually-averaged circulation are compared with the seasonallyvarying cyclo-stationary circulation in Section 7.2.

Offline time-dependent solver
To time-step the offline tracer-transport model we use 12 piece-wise constant monthly-averaged tracer-transport matrices.Within a given month the equation is stepped forward in time using an implicit first-order Euler-backward scheme for the vertical diffusion, and explicit schemes for the advection (3rd-order Adams-Bashforth), horizontal diffusion and source (first-order Euler forward): where i is a vector of ones and n n is a diagonal matrix formed from the elements of n n , which are the relative change in the surface cell thicknesses.The matrix-vector product Ai is the advective volumeflux convergence, and is non-zero only for the surface grid-cells.It does however average to zero over a full year because of the velocity corrections described in Section 4.1, ensuring that n is periodic with a period of one year.n 0 is initialized to g 0 =Dz 1 , where g 0 is the sea-surface height displacement saved from the parent model at the beginning of the one year run used to save the fields that are used to construct the transport operator./ n and / nþ1 are the tracer concentrations at time-step n and n þ 1. s is the source-sink function, which can depend on / n and n n .At the beginning of each month, the time-stepping scheme for the advective term is initialized using an Euler-forward step followed by a 2nd-order Adams-Bashforth step before switching to the 3rd-order Adams-Bashforth scheme.

Periodic equilibrium Newton-Krylov solver
To find cyclo-stationary equilibriums of the seasonally varying circulation model we apply Newton's method to GðxðtÞÞ ¼ 0, where GðxðtÞÞ where T ¼ 1 year.Given an initial iterate x 0 we obtain the improved iterate x 1 ¼ x 0 þ dx 0 by solving the following linear system where J @GðxÞ=@x is the Jacobian matrix.Unfortunately, J is a dense matrix because the integral in (13), evaluated with a timestep size of $3 h, allows a tracer perturbation isolated in one grid box at time t to propagate to every other grid box at time t þ T.
As a result, it is impossible to evaluate or even store the full Jacobian matrix in memory let alone invert it.To solve ( 14) we must therefore use a Krylov subspace method.For this we use a preconditioned generalized minimal residual method (with restarts) (GMRES) as implemented in the Matlab Newton-Krylov solver nsoli.mdescribed in Kelley (2003).

Preconditioner
To accelerate the convergence of the Krylov-subspace solver we follow Li and Primeau (2008) in using a preconditioner of the form where Q is a sparse matrix formulated in terms of the annual-average of the M matrix, defined in (4), As shown by Khatiwala (2008), the form of this preconditioner can be motivated as the inverse of a sparse approximation of the full Jacobian matrix J @GðxðtÞÞ @xðtÞ ; In the Newton-Krylov applications of Li and Primeau (2008) and Khatiwala (2008) the model resolutions were sufficiently coarse to make it possible to directly compute the LU factorization of Q .For the present model the factorization of Q requires nearly 256 GB of memory, making it impractical to use the LU factorization of the full matrix.Motivated by the desire to apply the solver to problems involving multiple coupled tracers and by the expectation that future versions of the model will have even higher resolution, we have developed a simple coarse-graining procedure for computing an approximate inverse of Q that requires much less memory.
The coarse-graining procedure consists of averaging successive 2 Â 2 horizontal blocks of grid-cells and then interpolating the resulting coarse-grained field back onto the full resolution grid by simply copying the tracer concentration in the coarse 2 Â 2 block into each of the four sub-gridboxes.This process of ''lumping'' tracer information onto a coarser grid and of ''spraying'' a lumped solution back to the full grid is illustrated in Fig. 3.It is achieved by constructing a pair of matrix operators, L and S defined as follows: where w denotes a column vector whose elements are the gridbox volumes, prime denotes matrix transpose, denotes the Kronecker product, . We refer to the N=4 Â N restriction matrix L as the lump operator because it lumps four neighboring grid-cell concentrations using a volume weighted average to produce a coarse grained concentration field.The N Â N=4 interpolation matrix S corresponds to the spray operator because it copies the concentration in the coarse grained grid box into each of the 4 sub-gridboxes of the full resolution model.After masking out the land points, the lump and spray operators are applied to the N Â N operator Q to reduce it to a $ N=4 Â N=4 matrix without appreciably changing its sparseness.The resulting preconditioner is then Note that the L and S matrices that pre-and post-multiply Q in (19) do not cancel with those that post-and pre-multiply Q À1 c in (20) because they are not square and therefore not invertible.
In summary, each Krylov iteration consists of time-stepping the offline-tracer transport model forward in time for one year, computing the drift in the tracer concentration between the beginning and the end of the one-year run, and then ''multiplying'' the resulting residual drift by P. In practice we do not explicitly form the matrix P or the inverse matrix Q À1 c .Instead we apply the preconditioner by lumping the resulting residual drift onto the coarser grid, applying the inverse of Q c by doing the backsolve with the precomputed lower and upper (LU) triangular factors, and spraying the result back onto the fine grid before subtracting the residual drift from the result.For tracers such as radiocarbon whose source terms depend only linearly on concentration, the LU factorization needs to be done only once.

Results
We applied the solver to the natural radiocarbon equations described in Section 2. We first describe the solver's convergence rate and the effect of seasonal variations in the transport operator on the annually averaged equilibrium D 14 C distribution.Then we present the resulting D 14 C distributions, which we compare to observations.

Convergence rate
Fig. 4 shows the convergence to the equilibrium solution expressed as the root-mean-square (rms) drift in the D 14 C for three cases: (1) The explicit time-stepping approach to equilibrium starting from an initial state obtained by interpolating the GLODAP gridded product onto the model grid.This case corresponds to the convergence rate that would have been obtained if radiocarbon had been simulated using the full parent model.(2) The approach to equilibrium using the implicit Newton-Krylov solver initialized using the interpolated GLODAP state and (3) the Newton-Krylov solver initialized using the steady-state solution based on the annually-averaged circulation.The latter approach would be particularly useful for applying the method to a tracer for which there are insufficient observations to use as an initial estimate.
Using time-stepping a root-mean-squared (rms) residual drift of 4 Â 10 À2 ‰=year was obtained after 100 model years of simulation.Using the Newton-Krylov solver initialized from GLODAP, the rms drift dropped to 4 Â 10 À4 ‰=year after 5 Newton iterations (corresponding to 23 years of simulation), at which point the OC-MIP-2 equilibrium requirements for D 14 C were met, i.e. more than 98% of the ocean volume had a D 14 C drift of less than 10 À3 ‰/year.After 7 Newton iterations (and a total of 66 years of simulation), the rms D 14 C drift was 7 Â 10 À12 ‰=year, and the Newton-Krylov solver stopped because it had reached our prescribed error tolerance of 10 À9 ‰=year.For the case where the Newton-Krylov solver was initialized using the steady-state solution, the OCMIP-2 equilibrium threshold was reached after 22 years of simulations.After 7 Newton iterations (and a total of 66 years of simulation), the rms D 14 C drift was 1 Â 10 À12 ‰=year.As expected, because of the linearity of Eq. ( 3), the rate of convergence was independent of the initial iterate used.
In terms of cpu time, the computational overhead associated with the Newton-Krylov solver (nsoli.m)was negligible compared to the cost of evaluating Eq. ( 13) i.e. of the cost of time-stepping the model forward in time for one year.In all our test cases we kept the GMRES restart parameter fixed at a value of 100, implying that the Newton-Krylov solver may store up to 100 three-dimensional tracer fields.This memory requirement, while substantial, is considerably less than the memory requirements needed to factor ($130 GB) and store ($77 GB) the LU factors of Q c used for the preconditioner.It is conceivable that further improvements in the efficiency of the solver could be achieved by increasing the restart parameter, but we did not explore this possibility.

L S ne mesh coarse mesh
Fig. 3.The Lump (L) and Spray (S) operators are used to implement the coarsegraining technique.The L operator transfers the tracer concentration on the fine grid to the coarse grid using a volume-weighted average of the fine-mesh gridboxes inside a coarse-mesh gridbox.The S operator transfers the tracer concentration on the coarse grid to the fine grid by copying the concentration in the coarse-mesh gridbox to each of the fine-mesh gridboxes inside the coarse-mesh gridbox.We also investigated other cases.If we initialize the explicit time-stepping method using the solution for the steady-state circulation model, the rms residual drift after 100 years is 1:2 Â 10 À2 ‰=year.The rate of convergence is similar to that shown for the case in which the model was initialized with the GLODAP state.Using the Newton-Krylov solver without a preconditioner, the rate of convergence is no better than that obtained from the time-stepping method.The rms of the residual after 100 years was 6:4 Â 10 À3 ‰=year, and there was no further progress during the subsequent 300 Krylov iterations, illustrating the critical importance of an effective preconditioner.
To demonstrate that the Newton-Krylov (NK) solution obtained using the offline model is an approximate equilibrium of the parent OGCM, we conducted two 20-year radiocarbon simulation in the parent model to compare the drift for the case where the model is initialized with the NK solution to the case where the model is initialized using the GLODAP pre-bomb radiocarbon.For these runs the dynamical state of the OGCM continued to evolve from the state at the end of the 150 year spin-up run.We used the same simplified representation for the air-sea exchange of 14 C that we used in the offline model.Fig. 5(a) shows the globally averaged rms D 14 C drift as a function of time for the two different initial conditions.For comparison we have reproduced on this plot the first 20-years of the rms D 14 C drift for the offline model initialized using the GLODAP data (see Fig. 4).The offline model and the parent OGCM show similar trends in the rms drift when initialized with GLODAP but there is nevertheless a significant difference that is most likely due to the approximation of the seasonal cycle in the offline model with a piecewise-constant transport operator with monthly resolution.The OGCM initialized with the offline model's equilibrium state has a substantial drift at the beginning of the run but decreases rapidly to less than 0.09‰/year at the end of the run.This drift is smaller than the GLODAP-initialized run by a factor of $10 at the beginning of the run and by a factor of $4 at the end of the 20 year run.The fact that the NK-equilibrium is only an approximate equilibrium to the parent model is most likely due to the poorly resolved seasonal cycle in the offline model.It should be possible to reduce this error by approximating the seasonal cycle in the offline model using two or more transport matrices per month.Another possibility would be to perform the one-year runs in the Newton-Krylov solver with the parent model instead of with the monthly transport matrices.In such a solver the transport matrix would be used only for the preconditioning step.In the deep ocean where the seasonal cycle is less important, we expect the offline model to better capture the parent model's circulation.Fig. 5(b) shows the averaged D 14 C value for the deep North Pacific (below 1000 m and north of 10°N) where the tracer equilibration times are longest.The upper panel shows the GLODAP-initialized run and the lower panel shows the NK-initialized run.The scale for the y-axis is the same for both panels except for an offset of À105‰.The run initialized with the offline model's equilibrium solution shows no discernable trend ( K 0:007‰/year) whereas the run initialized with GLODAP shows a downward trend of approximately 0:05‰/year.Extrapolating this trend gives an estimate of $2000 years for the model initialized using GLODAP to reach À329‰, the approximate equilibrium found by the offline NK solver.This estimate is likely to be an underestimate by a factor of two or more because the approach to tracer equilibrium is expected to follow an exponentially decaying function rather than a linear trend.
7.2.Importance of circulation seasonality on D 14 C Fig. 6 shows the differences in the D 14 C solutions obtained with and without the seasonality in the transport operator.The results are shown as a percent change of the average of the seasonallyvarying solution relative to the solution based on the annuallyaveraged transport operator.Positive percent changes indicate more depleted D 14 C values, i.e. older waters.With a few exceptions, discussed below, most of the differences are less than 10%.Overall, the bias from ignoring the seasonal variations in the circulation results in D 14 C values that are 1:5% less depleted.All basin zonal averages show increases in radiocarbon age (i.e. more depleted D 14 C relative to the solution based on the annually-averaged transport) in the upper part of the ocean between 100 m and 1000 m.The differences are proportionally largest in areas of strongest seasonal variability in the circulation.The difference between the solutions based on the cyclo-stationary and annuallyaveraged transport is especially notable in the northern part of the Indian Ocean where the overturning circulation changes sign on seasonal timescales (see Fig. 2).
In the Indian Ocean the difference between the cyclo-stationary solution and the steady-state solution penetrates all the way to the bottom of the basin with the cyclo-stationary solution having D 14 C values that are more than 30% more depleted in the depth range between 500 m and 1000 m.The surface of the Indian Ocean, on In the depth range between 100 m and 1000 m the Southern Ocean of the cyclo-stationary solution has D 14 C values that are 5-10% more depleted than the solution based on the steady time-averaged circulation.In the same depth range the North Pacific is also less well-ventilated in the cyclo-stationary solution compared to the solution based on the steady circulation with D 14 C values that are more depleted by at least 5%.estimate from GLODAP.Only the grid cells for which GLODAP data is available were used for the comparison.The light-gray band around the GLODAP values corresponds to an upper bound error estimate of AE30‰ (Gebbie and Huybers, 2012).A comprehensive error analysis of the GLODAP natural radiocarbon is not available, but sources of uncertainty for the GLODAP values are laboratorybased counting errors (Sabine et al., 2005); the separation of D 14 C values into natural and bomb-produced components (Rubin and Key, 2002), which includes regression errors of measured dissolved inorganic carbon (DIC) radiocarbon on potential alkalinity and on coral radiocarbon (Rubin and Key, 2002); and mapping errors (Key et al., 2004).The dark-gray band around the POP2 simulated values corresponds to the seasonal variation in the simulated D 14 C value.Interestingly, the error made in ignoring the seasonal variability in the circulation (see Fig. 6) is much larger than the amplitude of the seasonal cycle in the solution that properly accounts for the seasonally varying circulation.
In the Southern Ocean, south of 30°S, the simulated D 14 C is slightly more depleted than the observational estimate.The dip in the equatorial region in the GLODAP data is underestimated in the Atlantic and Indian basins, but not in the Pacific basin.In both the North Atlantic and the North Pacific, the simulated values are more depleted than the GLODAP values.Despite these differences, the simulated surface D 14 C values are mostly within the limit of uncertainty in the observations.In the North Pacific the model D 14 C is close to the uncertainty upper bound which indicates the model surface value is potentially in disagreement with the observations.
Our simulations used a very simple formulation for the air-sea exchange of radiocarbon.Further refinement of the air-sea gas exchange parameterization to include the effect of spatially varying air-sea CO 2 fluxes is expected to be important for the purpose of simulating the transient pulse of bomb radiocarbon (Krakauer et al., 2006), but for the purpose of simulating the prebomb radiocarbon component the simple parameterization appears adequate -the model's D 14 C values at the surface matched the observational estimates within their uncertainties in all the basins except perhaps in the North Pacific.The possible discrepancy in the North Pacific is likely due to the upwelling of deep waters in the model that are depleted by more than 80‰ compared to the observations and probably not the result of inadequacies of the simple air-sea gas-exchange parameterization.
Fig. 8 compares the simulated D 14 C values to the GLODAP estimates in the interior of the ocean.In agreement with the observations, the modeled D 14 C values are less depleted (younger) in the Atlantic, showing the substantial influence of ventilation from the Labrador Sea and the GIN Sea.In these North Atlantic waters the match between the model and GLODAP is excellent.However, the aging of waters as one moves southward in the model towards the Southern Ocean and northward into the Pacific and Indian Oceans is much greater in the model than in the GLODAP observations.In the North Pacific the difference between the model and GLODAP is most extreme, with the POP2 model D 14 C values reaching below À320‰ over a large portion of the water column, whereas the most-depleted GLODAP values are close to À240‰.This difference between the model and the observations is significantly greater than the 30‰ limit on the uncertainty in the GLO-DAP estimate.
Fig. 9 shows the model ventilation deficiencies using vertical D 14 C profiles separated by latitude band for each of the major ocean basins.Looking from the North Atlantic to the Southern Ocean (right to left in the figure), there is good agreement for water in the north end of the basin, where there is strong convection.There is a small more-depleted bias moving southward.This bias is strongest on the Eastern side of the basin as seen in Fig. 8.The sections of the Southern Ocean all have similar biases with the model D 14 C value becoming progressively overly-depleted with increasing depth.In the Indian Ocean the bias increases moving northward (left to right), with the largest biases at depth.Moving northward in the Pacific Ocean, the bias increases continuously, but unlike the case of the Indian Ocean, the maximum bias (an extremely depleted D 14 C value of À370‰) is centered above the bottom at a depth near 3000 m.The model's C14-age difference between the North Pacific at 3000 m and the surface of the Southern Ocean where the bulk of the North Pacific waters at 3000 m have last been in contact with the surface is $ 2600 years, which makes it more than two times bigger than the corresponding GLO-DAP value of only $ 1250 years.The bulge at 3000 m suggests that while the northward flow along the bottom is sluggish, southward flow in the mid-depth region of 3000 m further exacerbates the difference.

Discussion
We found that the model's North Atlantic D 14 C values matches well with the GLODAP data throughout the whole water column.The sill-overflow parameterization for the Denmark Strait and the Faroe Bank Channel, which were actively producing overflows in our simulation, appear to be successful at ventilating the deep North Atlantic.We speculate that the fact that the sill-overflow parameterization did not produce any overflows in either the Ross Sea or Weddell Sea is partly responsible for the D 14 C biases observed in the simulation -a more vigorous rate of bottom water formation in the Southern Ocean would help improve the ventilation of the Indian, Pacific and South Atlantic basins.It is not so clear what combination of changes to the surface forcing or to the subgridscale mixing parameterizations would be needed in order to produce a density profile that is favorable to initiating overflows in the Southern Ocean.Another likely cause for the weak contribution from the Southern Ocean to the model's deep water masses is likely due to the difficulty of representing in a one-degree resolution model the sea-ice interactions needed for deep water formation along the Antarctic continent.
It is interesting in this regard to compare the excessively depleted D 14 C values obtained with our model, À370‰, to the most-depleted D 14 C values obtained in previous modeling studies.Tables 1-3 list previous studies that presented D 14 C simulations done with ocean general circulation models.The tables also list the resolution of the models and whether or not the model used the GM parameterization.Several of the studies listed in the table ran multiple radiocarbon simulations with different forcing or different sub-gridscale parameterizations.Where more than one D 14 C value is given, it indicates the range of values that were found in the study; otherwise the authors' pick of the ''best'' simulation result is shown.
Table 1 lists studies that were made before the OCMIP-2 era.Most of the models are early versions of models that are in use today.All of them have horizontal resolutions that are coarser than 3 Â 3 and have between 10 and 21 levels in the vertical.Their most-depleted North Pacific D 14 C values range from À170‰ to À260‰.Contrary to our model results, the studies that did not match the observed estimate tended to produce deep North Pacific D 14 C minimums that were insufficiently depleted.The two models with the most depleted North Pacific D 14 C values, (Duffy et al., 1997;England and Rahmstorf, 1999), were produced with models that used the GM parameterization, which is also included in our version of POP2.Sensitivity experiments performed in these two studies showed that the GM isopycnal-layer thickness diffusion tended to limit the depth of deep convection in the model and as a result tended to slow the renewal rates of the deep Pacific Ocean waters.
Table 2 lists studies that were submitted to the OCMIP-2.See Doney et al. (2004), Orr (2002), Dutay et al. (2004);and Matsumoto et al. (2004) for studies focused on comparing the results from the OCMIP-2 models.The horizontal resolution for these models ranged from 5 to 2 in the zonal direction and from 0:8 to 5 in the meridional direction.The vertical resolution ranged from 15 to 30 levels.Excluding the PRINCE model, the most-depleted North Pacific D 14 C ranged from À190‰ to À250‰.The NCAR MOM1.1 model, a fore-runner of the POP model, had a maximum-depleted value of À230‰, one of the best matches with observations.
The PRINCE model D 14 C value of À340‰ was an outlier; however, the same modeling group also submitted an alternate version of the model with an improved eddy-diffusion parameterization that had a most-depleted North Pacific D 14 C value of À220‰, which is in reasonable agreement with the observational estimate.
This alternate version of the model described in Gnanadesikan et al. (2004), showed that a proper calibration of the model's eddy-diffusive parameterizations could produce acceptable D 14 C distributions without negatively impacting the simulation of other tracers such as temperature, salinity, oxygen, phosphate and CFC-11.Table 3 lists more recent studies performed after the OCMIP-2 era.These models have a broad range of resolution ranging from 10 Â 10 to 1:4 Â 1:4 in the horizontal, and with vertical resolutions ranging from 7 to 32 levels.All but one of these studies produced minimum North Pacific D 14 C values that are in agreement with the GLODAP estimate.Interestingly, the one model that produced a minimum D 14 C value that is outside the observational error estimate, (Graven et al., 2012) used the same POP2 model as our study, but configured with a substantially coarser vertical and horizontal grid.
Taken together, these studies show that coarse resolution models were generally able to simulate prebomb D 14 C distributions that agreed reasonably well with the observationally based estimate.We believe that higher resolution models should be able to do as well or better, provided that modelers have the ability to perform enough radiocarbon simulations to properly calibrate sub-gridscale parameterizations.Doing so should lead to improvements in the representation of deep-ocean ventilation processes in their models.

Conclusions
We have presented an offline implicit solver for obtaining the equilibrium prebomb radiocarbon distribution in an ocean circulation model of moderately-high resolution.The solver uses a Newton-Krylov method with a preconditioner based on a direct LU factorization of a coarse-grained annually-averaged transport matrix.We expect that the coarse-graining method we presented here will provide a general strategy for developing effective preconditioners that will make it possible to apply computationallyefficient Newton-Krylov solvers to models of increasing resolution.While tuning the physical parameters of an ocean model will require on-line spin-up of the dynamics of the parent model, the implicit solver should be useful for providing ventilation diagnostics of the deep ocean that can be computed without the need for long runs of the parent dynamical model.We have shown that this solver makes it practical to use D 14 C as a diagnostic for evaluating the deep-ocean ventilation of ocean models currently being developed for coupled ocean-atmosphere climate models.
For our application of the solver to the CESM POP2 model we found that the simulated D 14 C values agree well with the GLODAP estimate at the sea surface and in the North Atlantic, but that there are significant biases in the deep North Pacific with D 14 C values that are overly depleted by up to $110‰.This bias is larger than that obtained in all previously published radiocarbon simulations even though the resolution of this model is better than that used in any previously published prebomb radiocarbon simulation.
While we believe that increasing the resolution of ocean models is necessary for improving the circulation of ocean models, our results show that improvements in fidelity with increasing resolution are not automatic.Substantial efforts must be dedicated to calibrating sub-gridscale eddy-diffusive parameterization and to improving surface forcing fields.Natural radiocarbon provides a useful metric with which to evaluate the ventilation of the deep ocean in models and our implicit solver makes it computationally feasible to apply this metric.

Fig. 2 .
Fig. 2. Annual, DJFM, and JJAS averages for the global (left) and for the Atlantic (right) meridional overturning circulation (MOC), in Sverdrup (10 6 m 3 s À1 ).Red colors indicate a clockwise overturning.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig.4.Convergence rate for D 14 C approach to equilibrium for the cases of: (1) a usual time-dependent simulation initialized from the GLODAP observations interpolated onto the model grid (red-squares); (2) using the Newton-Krylov solver with the same initial iterate as used for the time-dependent case (black asterisks); and (3) using the Newton-Krylov solver with an initial iterate obtained by directly solving the steady-state of the model based on the annually averaged transportoperator (blue pluses).After 5 Newton iterations (equivalent to 22 model years in case 3, and 23 model years in case 2) the OCMIP-2 equilibrium criteria has already been met.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .
Fig. 5. (a) Globally averaged rms drift for D 14 C in a 20-year run of (i) the parent OGCM initialized using the offline equilibrium solution obtained with the NK solver (blue þ), (ii) of the parent OGCM initialized using the GLODAP prebomb radiocarbon interpolated on the model grid (black squares), and (iii) of the offline model also initialized with GLODAP (red squares).(b) The D 14 C in the deep North Pacific below 1000 m and north of 10°N as a function of time for the OGCM runs initialized with GLODAP (upper panel) and with the offline equilibrium solution (lower panel).Note that the y-axis for the upper and lower panels have the same scale but are offset by À105‰.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. D 14 C in the surface layer versus latitude.The black solid line is the GLODAP data; the dashed line is the annually-averaged cyclo-stationary solution.The light gray bands are AE30‰ from the GLODAP data.The dark gray band around the dashed line represents the range in the seasonal cycle values.The panels show the average by latitude for the global ocean and for each major basin.

Fig. 9 .
Fig. 9. Depth profiles of D 14 C in (‰) by latitude bands, in the Atlantic, Indian, and Pacific.The solid black line is GLODAP data; the dashed line is the simulated value.The light gray band is AE30‰ from the GLODAP value.The dark gray band around the dashed line gives the seasonal variation in the value: this is nearly invisible on this scale.

Table 1
Ocean model studies before OCMIP-2 whereD 14 C was used to validate or study the model.''GM?''states whether or not the Gent-McWilliams parameterization was used.

Table 2
Orr (2002)ral circulation model studies participating in OCMIP-2 whereD 14 C results were submitted for the comparison.The references give descriptions of the models.If the reference does not contain the D 14 C results, then the result for the most-depleted Pacific D 14 C value was obtained fromOrr (2002)summary of the results.If more than one value is given as the most-depleted Pacific D 14 C value, the study produced the range of results given.Model with a terrain-following coordinate system is marked by y.

Table 3
Ocean model studies post-OCMIP-2 where D 14 C was used to validate or study the model.Isopycnal models are marked with a y.