On the time to tracer equilibrium in the global ocean

. An important issue for the interpretation of data from deep-sea cores is the time for tracers to be transported from the sea surface to the deep ocean. Global ocean circulation models can help shed light on the timescales over which a tracer comes to equilibrium in different regions of the ocean. In this note, we discuss how the most slowly decaying eigenmode of a model can be used to obtain a relevant timescale for a tracer that enters through the sea surface to become well mixed in the ocean interior. We show how this timescale depends critically on the choice between a Neumann surface boundary condition in which the ﬂux of tracer is prescribed, a Robin surface boundary condition in which a combination of the ﬂux and tracer concentration is prescribed or a Dirichlet surface boundary condition in which the concentration is prescribed. Explicit calculations with a 3-box model and a three-dimensional ocean circulation model show that the Dirichlet boundary condition when applied to only part of the surface ocean greatly overestimate the time needed to reach equilibrium. As a result regional-“injection” calculations which prescribe the surface concentration instead of the surface ﬂux are not relevant for interpreting the regional disequilibrium between the Atlantic and Paciﬁc found in paleo-tracer records from deep-sea cores. For


Introduction
In a recent article Wunsch and Heimbach (2008) (hereinafter referred to as WH08) pose an important question for the interpretation of paleoceanographic proxy tracer data: "How long to oceanic tracer and proxy equilibrium?" One specific motivation for posing the question is the large difference in the apparent time for the δ 18 O of deep Atlantic and Pacific oceans to reach equilibrium subsequent to the input of depleted δ 18 O water from the melting of ice sheets noted by Skinner and Shackleton (2005). WH08 suggest that a transient tracer disequilibrium can more than account for the 3900 year lag in the Pacific relative to the Atlantic recorded in the sediment. Indeed, simulated transient tracer experiments performed by WH08 suggest that in the case where a tracer enters the ocean over a limited-area patch of the high latitude Northern Atlantic or Southern Ocean, the deep Pacific can lag the deep Atlantic by as much as 4000-6000 years before it reaches 90% of its equilibrium value. While we agree entirely with the main point of the WH08 articlethat a substantial passage of time is needed before a tracer becomes uniformly mixed throughout the ocean after it enters through the sea surface, and that this transient response can result in significant regional differences in the time to reach uniformity -we wish to bring attention to the importance of the choice of boundary condition on the implied time to reach equilibrium.
The long equilibration timescale obtained by WH08 depends critically on their choice of surface boundary condition. In their numerical simulations they prescribe the tracer Published by Copernicus Publications on behalf of the European Geosciences Union.  The injection is uniformly distributed in a 50 m thick patch with an area 1.4×10 7 km 2 located in the high latitude Northern Atlantic as indicated by black region in the inset map. The dashed curve corresponds to the Dirichlet boundary condition in which the tracer concentration is prescribed to unity for t>0 in the same 50 m thick patch used for the Neumann problem. At time t=0 the concentration of both tracers is zero everywhere in the ocean.
concentration in a patch at the surface of the ocean. If instead of prescribing the concentration of tracer they had prescribed the flux of tracer, the time to reach equilibrium would have been drastically shorter. This is illustrated in Fig. 1 which shows the time evolution of the concentration in the deep Atlantic ocean at (9 • N, 19 • W, depth=3195 m) of two tracers simulated using a three-dimensional ocean circulation model to be described in Sect. 3. One tracer, denoted by the dashed red curve, was introduced using the same type of boundary condition used by WH08. Its concentration was prescribed to be unity starting at t=0 over a 50 m thick patch in the high latitude North Atlantic Ocean. As t→∞ and all the ocean's water parcels have cycled through the surface patch, the concentration of this tracer approaches unity everywhere in the ocean. The other tracer, denoted by the solid blue curve, was introduced by prescribing an instantaneous uniform flux at time t=0 into the same 50 m thick patch. For this second tracer, the amount of tracer flux into the ocean was prescribed to be exactly the amount needed such that once the tracer is mixed uniformly throughout the ocean its concentration is equal to unity everywhere. By construction, both tracers have the same initial condition at time t=0 and the same equilibrium condition as t→∞. Their approach to equilibrium, however, is dramatically different. After 25 000 years of simulation the tracer with a prescribed concentration boundary condition has barely reached 90% of its final equilibrium concentration but the tracer that is introduced using a prescribed flux boundary condition is already in approximate equilibrium after only 1200 years.
The reason for the vastly different equilibration timescales is easy to understand physically. For the case in which the concentration is prescribed, the tracer is introduced into the ocean only gradually because the flux of tracer through the surface patch depends on the tracer gradient normal to the surface. As the tracer enters the ocean its gradient near the surface, and hence its flux through the surface, decreases rapidly because the tracer naturally spreads first to waters that are near the patch. Further uptake of tracer depends on the rate at which pristine waters can cycle through the patch so as to maintain a concentration gradient normal to the patch. For the case in which the flux of tracer is prescribed, all the tracer is introduced into the ocean instantaneously at time t=0 and the role of the circulation is then to mix the tracer throughout the ocean. Equilibrium is reached when the tracer becomes uniformly mixed everywhere. In case the reader objects that in reality tracers enter the ocean more gradually, we demonstrate in Sect. 3.2.1 that the overshoot seen in Fig. 1 for the case in which the flux is applied impulsively is greatly diminished if the flux is distributed over a thousand years or more. However, an important point for the interpretation of paleo-proxy records is that time lag for different points of the deep ocean to reach equilibrium when the flux is distributed in time is still controlled by the equilibration time of the impulsive response and remains drastically shorter than for the prescribed concentration response.
The boundary condition used by WH08 can be thought of as labeling fluid elements with a tracer as they circulate through the surface patch. The concentration, C, increases monotonically from C=0 at time t=0 to C=1 as t→∞. The tracer concentration, C(r, t), can therefore be interpreted as the fraction (by volume) of the water parcel centered at r at time t that has circulated through the patch at least once since the time t=0. The time t 90 (r) at which C=0.9 for the first time and used by WH08 as a measure of tracer equilibration time corresponds to the time at which 90% of the fluid elements at r have made contact with the surface patch since t=0. To the extent that the circulation can be thought of as being stationary, C(r, t), can also be interpreted as a cumulative distribution function of times, t, since the fluid elements at r have last circulated through the surface patch. Note that when the surface patch covers the entire sea surface the distribution is sometimes called the age distribution (e.g. Primeau, 2005). With this interpretation, t 90 (r), is the 90th percentile of the distribution function of times since fluid elements at r have made last contact with the patch. As is evident from Fig. 1, the time for an initially concentrated patch of tracer -such as would occur from a localized input of anomalous δ 18 O water from melting ice sheets -to mix F. Primeau and E. Deleersnijder: Tracer equilibrium 15 uniformly throughout the ocean is much shorter than the time needed for 90% of the fluid elements to cycle through the input patch. Consequently, the t 90 criterion used by WH08 greatly overestimates the time needed for a pulse of anomalous δ 18 O water to reach equilibrium.
Our goal in the following sections is to further elucidate the dependence of the equilibration time on the choice of boundary condition. However, one of the main conclusions from this article should already be evident -the relevant boundary condition for studying the differences in equilibration time between the Atlantic and Pacific ocean in the δ 18 O record is not the one used by WH08 in which the concentration is prescribed, but one in which the flux of tracer is prescribed. As we will discuss later, a prescribed concentration boundary condition can be used to extract relevant timescales for tracers that enter the ocean through air-sea exchange provided the air-sea gas exchange rate is sufficiently fast, but the prescribed concentration boundary condition must be applied over the entire air-sea interface and not as was done in WH08 over a patch of limited area with no-flux boundary conditions applying over the rest of the ocean.
A natural timescale for characterising the rate at which a system approaches equilibrium is the e-folding decay timescale of the most slowly decaying eigenmode of the system. As t→∞, the most slowly decaying eigenmode dominates the time evolution of the system because the other eigenmodes eventually decay sufficiently to make a negligible contribution. The spatial pattern associated with the most slowly decaying eigenmode together with the mode's decay rate provide a concise description of the regional differences in the equilibration time. The plan for this paper is therefore to compute the eigenmodes of a three-dimensional OGCM, and explore the dependence of the characteristic equilibration time on the size of the patch in which the tracer is introduced. Before presenting the results for the OGCM, we first illustrate in the simplest possible terms the main difference for the equilibration time of a problem with a prescribed flux (Neumann boundary condition) and a prescribed concentration (Dirichlet boundary condition) using a simple 3-box model. The available analytical solution for the simple 3-box model will allow the interested reader to verify qualitatively the results which we obtain numerically for the full OGCM.

A simple 3-box model example
To illustrate in the simplest possible terms how the time to reach equilibrium depends critically on the choice between a Dirichlet and Neumann boundary condition, we now consider the evolution of a passive tracer in a simple 3-box model. A detailed formulation of the governing equation for the 3-box model shown in Fig. 2 is presented in Appendix A. Expressed in matrix form, the governing equation is . 2. (a) Geometry of the 3-box ocean model where V 1 , V 2 and V 3 are the box volumes, αS and (1−α)S with 0≤α≤1 are the areas of the interfaces separating box 1 from 3 and box 2 from 3 respectively, and where S is the area of the ocean-atmosphere interface.
(b) Diagram showing the notation used for the tracer concentrations c 1 , c 2 and c 3 , for the fluxes between the boxes φ 1,3 =(c 1 −C 3 )αvS, and φ 2,3 =(c 2 −c 3 )(1−α)vS, and for the flux from the atmosphere or runoff, φ a,1 . where, is the model's flux-divergence operator (also known as the transport operator), where the state vector c=[c 1 c 2 c 3 ] represents the tracer concentration in each box and where s=[φ a,1 0 0] is a source term. The prime denotes vector transpose. The parameters in A that are not defined in the caption for Fig. 2 are the inverse timescale γ ≡vS/(V 1 +V 2 ), and the small parameter ≡(V 1 +V 2 )/V 3 .

Neumann boundary condition
We first consider the case in which a pulse of tracer is injected instantaneously into box "1" at time t=0 prior to which the tracer concentration was zero everywhere. The total amount of tracer injected is such that as t→∞ and the mixing of the tracer is complete, the equilibrium tracer distribution is one in which c 1 =c 2 =c 3 =1. The total amount of tracer injected at time t=0 must therefore be V 1 +V 2 +V 3 , and the concentration in box 1 at time t=0 + (immediately after tracer injection) is For t>0 there is no source or sink of tracer into the ocean so that the total amount of tracer is conserved. The time evolution of the tracer can then be obtained by propagating the initial tracer distribution forward in time with the exponential of the matrix A (e.g. Hirsch and Smale, 1974): where c o ≡c(t=0 + )=[(1+ )/(α ), 0, 0] is the tracer distribution immediately after injection when all the tracer is concentrated in box "1". Expanding the exponential in terms of the eigenmodes of A, we obtain the solution The first mode corresponds to a zero eigenvalue and is therefore independent of t. It captures the equilibrium state of the system which corresponds to a state in which the pulse of tracer is mixed uniformly throughout the three boxes. The last two modes decay exponentially with e-folding timescales given by the reciprocal of the eigenvalues, λ + = −γ and λ − = −(1 + )γ .
These two decaying modes describe the transient state of the system. The longest lived of the two modes with an e-folding timescale 1/λ + controls the length of time for the system to come to equilibrium. The adjustment timescale for the Neumann problem is therefore It is important to note that unlike the adjustment timescale τ D obtained for the Dirichlet problem to be considered in Sect. 2.2, τ N is independent of the relative sizes of the two surface boxes. As time progresses mixing and diffusion erases the details of the tracer injection so that the size of the initial injection box becomes irrelevant.

Dirichlet boundary condition
The box model analogue of the Dirichlet boundary condition considered by Wunsch and Heimbach (2008) consists of holding the tracer concentration in box "1" to unity starting at t=0. For t<0 the tracer concentration in each box is prescribed to be zero. Because c 1 is prescribed, the differential equation corresponding to the first row of A in Eq. (1) can be ignored and we are left with only the last two equations. The resulting governing equation for the Dirichlet boundary condition reduces to where c o ≡[c 2 (t), c 3 (t)] , and The solution to this system can be written as the sum, c o =c f +c h , of a steady forced solution and a transient homogeneous solution which can be written in terms of the exponential of a matrix Note that at time t=0, c h =−c f so that c o satisfies the zero initial condition. Expanding the matrix exponential into the eigenmodes of A oo , the time evolution of the tracer subject to the Dirichlet boundary condition can be expressed as follows where Both λ + and λ − are real and negative so that the homogeneous solution decays to zero and the equilibrium solution corresponds to the forced solution which is a state in which the tracer concentration is uniform throughout the boxes. The approach to equilibrium for the Dirichlet problem considered here is governed by the eigenmodes of the matrix A oo . The time to reach equilibrium is given by the mode whose e-folding timescale is longest, i.e.
Note that τ D depends on α, the relative size of the surface boxes. This is consistent with the results of WH08 who also noted that the timescale to reach equilibrium increases as the size of the surface patch where the tracer enters the ocean decreases.

Discussion of the simple box model results
If we form the ratio of τ D to τ N we obtain, which shows clearly that the timescale for the Dirichlet problem is much greater than the adjustment timescale for the corresponding Neumann problem.
The Dirichlet boundary condition implies a definite timedependent "air-sea" flux of tracer into box "1". Indeed, a consideration of the total tracer budget yields φ a,1 (t)=φ 1,3 =αvS [1−c 3 (t)]. From which we see that as c 3 (t) approaches its equilibrium value the implied flux of tracer decreases to zero. We see also that the implied flux is proportional to αS the area of the patch where the tracer is introduced. The smaller the area where the tracer is introduced into the ocean the more slowly the tracer is fluxed into the ocean. The slow equilibration timescale for the Dirichlet in comparison to the Neumann problem can be understood in terms of the slow rate at which the tracer is introduced into the ocean in the Dirichlet solution. This slow input rate is not relevant for interpreting the δ 18 O record because the timing and amplitude of the melt water pulse is causally independent from the concentration of δ 18 O-depleted water in the ocean. The melting of the ice prescribes a flux of δ 18 O-depleted water into the ocean and, given the melting rate, does not depend on the circulation as is the case for the Dirichlet boundary condition. The surface concentration of δ 18 O-depleted water on the other hand is determined by the balance between the prescribed flux and the mixing of surface and interior waters. It is therefore the relatively fast timescale of the Neumann problem that is the relevant one for interpreting the δ 18 O record.
An alternative way of thinking about the slow adjustment timescale of Dirichlet problem is to think of the solution to the Dirichlet problem as the concentration of surface water (e.g. Deleersnijder et al., 2002), i.e. the concentration of water particles that have touched box-"1" at least once. The smaller we make α the less likely it is for a water particle to have hit the box-"1" air-sea interface and thus have been transformed into a "surface water particle". The t 90 timescale defined by WH08 as the time at which the tracer concentration at some point r in the interior of the ocean equals 0.9 is precisely the time at which 90% of the water in a parcel at r has been in contact with the patch where the tracer is introduced. The solution to the Neumann problem shows clearly that a pulse of tracer becomes well mixed long before 90% of the water particles have made contact with the patch where the tracer is injected.
The box model is sufficient to demonstrate clearly how the Dirichlet versus Neumann boundary conditions lead to highly different adjustment timescales. However, in order to get more quantitative results relevant to the real ocean and to look at the relative disequilibrium between the Atlantic and Pacific Oceans, we now turn to a three-dimensional ocean circulation model.

Eigenmode analysis for a 3-D ocean tracer transport model
In order to compare the eigenmodes of a three-dimensional global ocean circulation model with a Neumann bound-ary condition (prescribed flux) to a model with a Dirichlet boundary condition (prescribed concentration), we use the three-dimensional tracer transport model of Primeau (2005).
The advantage of this model is that its advection-diffusion tracer transport operator is available in matrix form making it possible to compute the model's three-dimensional eigenmodes by solving a matrix eigenvalue problem. The tracer transport model is driven by the velocity field and eddy diffusion tensor field derived from a dynamical ocean general circulation model (OGCM  McWilliams, 1990) isopycnal eddymixing scheme. Convection is parameterized by using an implicit vertical diffusion scheme with a large diffusivity coefficient. Second-order centered differences are used on a ∼3.75 • ×3.75 • grid with 29 levels ranging in thickness from 50 m near the surface to 300 m near the bottom. The model has a velocity field and transport characteristics typical of OGCMs with similar resolution and produces a maximum Atlantic meridional overturning streamfunction of 18 Sv. It has been used to simulate transient and equilibrium tracers such as 14 C and CFC's (Krakauer et al., 2006) and biological tracers such as phosphate, dissolved inorganic carbon, and alkalinity (Kwon and Primeau, 2006;Kwon and Primeau, 2008). The water mass ventilation properties of the annually averaged circulation are described in detail in Primeau (2005), Primeau and Holzer (2006), Primeau (2006, 2008). The tracer transport model uses a state vector c of dimension 63 090×1 whose elements correspond to the tracer concentration in each ocean grid box to represent a threedimensional tracer concentration field. The advectiondiffusion flux-divergence operator subject to no-flux conditions on all boundaries, once discretized, can be expressed as 63090×63090 sparse matrix, A=∇·[u−K·∇], in which u is the fluid velocity and K is the eddy-diffusivity tensor.

Dirichlet boundary condition: prescribed tracer concentration
We present the problem formulation in discrete form. The interested reader is referred to Appendix B for the continuous formulation of the eigenmode problem. For the Dirichlet boundary condition the tracer concentration is prescribed in the patch so that only the equations for the grid points outside the patch need to be solved for. The in-patch points for which the concentration is prescribed to be unity appear on the right hand side as a forcing term to a reduced system of equations for the out-of-patch grid boxes, F. Primeau and E. Deleersnijder: Tracer equilibrium in which 1 s is a vector of ones of length N s , the number of grid points inside the patch. In writing Eq. (16), we used the partitioning of the state vector and matrix transport operator into in which s indicates the set of indices corresponding to model grid boxes inside the surface patch and o indicates the set of indices corresponds to grid boxes outside the patch. The eigenmode problem for the discretized model with a Dirichlet boundary condition consists of looking for solutions of the form c o (t)=ve −λt . Substituting this modal form into the homogeneous counterpart of Eq. (16) leads to the following matrix eigenvector problem, where N o is the number of grid boxes outside the patch. In the following we will assume that the eigenmodes are ordered in terms of their real part {λ n }≡τ −1 n such that τ 0 >τ 1 ≥τ 2 ≥ · · · ≥τ N .
Note that by eliminating the rows and columns corresponding to grid points inside the patch we have made the matrix in Eq. (18) depend on the location and size of the patch. This should be contrasted with the Neumann eigenmode problem. As we will see in Sect. 3.2, (Eq. 23) the eigenmode problem for the Neumann problem uses the full matrix and is therefore completely independent of the patch size and location. The e-folding decay rate of the most slowly decaying eigenmode for the Dirichlet problem will therefore depend on the size and location of the patch. This is consistent with the results of WH08 who found that the t 90 timescale tended to increase as the patch size was made smaller.

Approach to equilibrium for the Dirichlet problem
The solution of Eq. (16) can be expressed as the sum of a forced and a homogeneous solution. The forced solution describes the equilibrium state of the system and the homogeneous solution describes the transient approach to equilibrium. If we expand the homogeneous solution in terms of the system's eigenmodes, the full solution can be written as follows where the leading 1 o corresponds to the forced response due to the prescribed concentration of unity in the patch. We have separated out the leading eigenfunction, v 0 e −t/τ 0 , because its eigenvalue is purely real. To see why this is so, note that in the limit of t→∞ the time evolution of any tracer distribution that has an initial projection on v 0 will eventually become dominated by v 0 because all the other eigenmodes decay to zero more quickly. In this asymptotic regime, the tracer concentration would change sign at certain values of t if ω 0 =0. Because such changes of sign are a physical impossibility, ω 0 must necessarily vanish. We can also conclude from this argument that the leading eigenvector, v 0 is necessarily sign definite otherwise an initially positive tracer distribution would ultimately produce negative tracer concentrations once the other eigenmodes had decayed away. Furthermore, any positive tracer distribution must have a projection onto v 0 because in a connected domain with finite diffusivity all the other eigenmodes, v n for n=1, 2, · · · ∞ cannot be sign definite, that is they must have both positive and negative regions. This fact follows from the bi-orthogonality property of the eigenvectors of A and its adjoint which are both proper advection-diffusion tracer transport operators. The fact that any positive tracer distribution must have a projection onto v 0 for the Dirichlet problem means that v 0 e −t/τ 0 is always the relevant eigenmode for describing the asymptotic approach to equilibrium of the Dirichlet problem. Because the type of boundary condition that is used inside the patch (prescribed concentration) is different from the type of boundary condition that is used outside the patch (no flux) there is no single eigenspectrum relevant to each Dirichlet problem corresponding to a different patch size and location. A separate eigenvalue problem must be solved every time the surface patch is changed. We note also that because different types of boundary conditions are used inside and outside the patch the principle of superposition does not apply if one wants to add up the concentrations from different patches. This points again to the fact that the regional-"injection" runs of WHO8 cannot be used in a simple way to interpret real tracers for which we expect to be able to add the contribution from different patches to obtain the total tracer concentration.
To explore the dependence of the e-folding decay timescale of the most slowly decaying mode on the location and areal extent of the patch where the tracer concentration is prescribed, we randomly chose 20 points at the surface of the ocean and for each point we constructed a set of 6 patches ranging in area from that of a single grid box to that of the full ocean surface. In this way a total of 101 distinct surface patches were constructed. For each of these patches we constructed a matrix A oo by deleting the rows and columns of A corresponding to grid points inside the patch. We then solved the matrix eigenvalue problem in Eq. (18) for each of the 101 A oo matrices using the sparse matrix eigenvalue solver Arpack (Lehoucq and Sorensen, 1996) as implemented in Matlab's eigs function. For each case we recorded the efolding decay timescales, τ 0 =1/λ 0 and τ 1 =1/ {λ 1 }, of the two most slowly decaying eigenmodes. Figure 3 shows a scatter plot of τ 0 and as a function of the reciprocal of the area A of the surface patch. The figure shows clearly how the equilibration timescale for the Dirichlet problem increases without bound as A approaches zero. There is also a conspicuous scattering of τ 0 for small area patches indicating a large sensitivity to the precise location of the patch. This sensitivity is expected due to local differences in the vertical velocity, diffusive mixing and penetration depth of the mixed layer. In contrast, the variability in timescale τ 1 (not shown) associated with the decay rate of the second most slowly decaying eigenmode decreases and approaches the decay rate of the most slowly decaying eigenmode of the Neumann problem as A→0. This behavior of τ 0 and τ 1 is consistent with the fact that as A→0 the eigenspectrum of the Dirichlet problem approaches the one for the Neumann problem. In other words, as A→0 the number of grid boxes inside the patch decreases until A oo →A and 1/τ 0 →0 as is required for the Neumann problem where the no-flux boundary condition ensures that the total amount of tracer is conserved.

Neumann boundary condition: prescribed flux
We again present the problem formulation in discrete form.
The interested reader can refer to Appendix C for the continuous formulation of the eigenvalue problem with Neumann boundary conditions. The tracer transport problem in which a flux of tracer is prescribed through a given surface patch at time t=0 can be written in terms of a matrix system of differential equations of the form subject to the initial condition where s is a vector with ones in the elements corresponding to grid boxes inside the patch and zeros in the elements corresponding to grid boxes outside the patch and where δ(t) is the Dirac-delta function. The dimensionless scalar ρ ensures that enough tracer is injected in the patch at time t=0 so that the asymptotic equilibrium tracer concentration is equal to unity everywhere. The parameter ρ is therefore equal to the ratio of the total volume of the ocean to the volume inside the patch. Equation (20) can be reduced to a homogeneous initial value problem (i.e. with no source term on the right hand side) by integrating the equation from t=− to t= and letting →0, to obtain d dt c + Ac = 0, in which the effect of the pulse of tracer injected into the ocean at time t=0 is encoded in the problem's initial condition. The eigenmode problem for the discretized model with Neumann boundary condition consists of looking for solutions of the homogeneous equation that are of the form c(t)=ve −λt . Substituting the modal form into the homogeneous equation leads to the following matrix eigenvector problem, Av n = λ n v n , for n = 0, · · · , N − 1, where N =63 090 is the number of grid boxes in the model. The eigenvectors, v n , capture the spatial pattern of the eigenmodes and the eigenvalues, λ n are such that the reciprocal of their real parts give the e-folding decay rate of the corresponding eigenmodes. In the following we will assume that the eigenmodes are ordered in terms of their real part {λ n }≡τ −1 n such that τ 0 >τ 1 ≥τ 2 ≥ · · · ≥τ N .

Approach to equilibrium for the Neumann problem
The solution to Eq. (20) is obtained by projecting the initial condition, c=ρs onto the eigenmodes of A, c(t) = e At (ρs) = N −1 n=0 a n v n e −λ n t = 1 + N −1 n=1 a n v n e (−1/τ n +iω n )t , where N is again the total number of grid boxes and hence eigenmodes, and where the leading vector of ones in the second line corresponds to the constant eigenfunction with the  zero eigenvalue. All the eigenmodes with n≥1 have a positive τ n and therefore decay exponentially with time. As t→∞ only the constant eigenmode with a zero eigenvalue survives to produce the asymptotic equilibrium state. In general, the approach to the asymptotic tracer field as t→∞ is dominated by v 1 , because the other eigenmodes decay more quickly. In the unlikely situation where the initial tracer distribution does not project onto v 1 , the asymptotic tracer field will be dominated by v 2 e −λ 2 t , the second most slowly decaying eigenmode and the approach to equilibrium will be even faster.
To compute the slowly decaying part of the eigenspectrum of A we again used Arpack (Lehoucq and Sorensen, 1996) as implemented in Matlab's eigs function. As required by tracer conservation, A has a constant eigenvector with a zero eigenvalue, λ 0 =0. The most slowly decaying eigenmode for our model has eigenvalue, λ 1 ≡1/τ 1 −iω 1 , which corresponds to an over damped mode with an e-folding decay time of τ 1 =629 years and a period of T n =2π/ω 1 =5190 years. The next most slowly decaying eigenmode has an e-folding decay timescale of τ 2 =365 years. Figure 4 shows a contour plot of the phase and amplitude of the corresponding eigenfunction, v 1 in the deep ocean.
For our OGCM, the approach to the uniform tracer distribution is governed by an over-damped exponentially decaying mode, v 1 e (−1/τ 1 +iω 1 )t . By factoring out e −t/τ 1 in the solution's eigenmode expansion in Eq. (24), we can rewrite the solution as follows c(t) = 1+e −t/τ 1 a 1 v 1 e iω 1 t + N −1 n=2 a n v n e (−1/τ n +1/τ 1 +iω n )t , which shows that the exponentially decaying regime is itself approached with a timescale determined by the relative decay rate of the next most slowly decaying eigenmode through the formula τ r ∼τ 1 τ 2 /(τ 1 −τ 2 ). For our model is τ r =870 years.
Assuming that the initial projection on v 1 is not zero, as is the case for our model, the relative disequilibrium between points at two different locations in the deep ocean, in grid boxes i and j say, can be obtained for times t τ r from the spatial structure of v 1 alone. If we consider two points in the ocean at grid boxes i and j then the asymptotic disequilibrium of these two points is where v 1i and v 1j are the ith and j th elements of v 1 . The time lag t (i, j ) between the equilibration time of boxes with indices i and j can be estimated from the requirement that which leads to the following formula, provided neither v 1i or v 1j vanishes. At the few points where either v 1i or v 1j vanishes, the next eigenmode in the expansion (which decays even more rapidly) would need to be taken into account. For our model, the lag in the equilibration time between any two points at depths below 3100 m is always less than 1200 years and much less on average. These time lags are much smaller than those obtained from the boundary condition used by WH08, i.e. Dirichlet boundary condition on regional patches. Figure 5 contrasts the time evolving response at two points in deep ocean to a Dirac-delta function pulse of tracer in the surface ocean. One point is in the Atlantic and the other point is in the Pacific. The tracer concentration crosses its equilibrium value before 2000 years at both sites and is pretty much in equilibrium at about 3000 years. The Pacific lag with respect to the Atlantic is generally less than 1000 years. Figure 6 shows the response at the same two points for the case where the injection of tracer is distributed in time according to a Gaussian pulse with a standard deviation of 2000 years. Distributing the tracer flux in time as opposed to injecting it instantaneously using a Dirac-delta function has the effect of reducing the overshoot past the asymptotic equilibrium value and tends to further decrease the lag between the Pacific and Atlantic responses.

Robin boundary condition: prescribed linear combination of flux and concentration
In the case where the tracer enters the ocean through air-sea gas exchange it is often the case that the information we have is for the atmospheric concentration of the tracer. The appropriate boundary condition is then one in which the flux of tracer through the air-sea interface is taken to be proportional to the difference in the concentration of the tracer between the bulk ocean and the bulk atmosphere (e.g. Krakauer et al., 2006). This implies a surface boundary condition of the form, where F is the gas flux out of the ocean; C s is the concentration of the gas in surface water; C a is the surface ocean concentration of the gas in equilibrium with the partial pressure p of the gas in the air over the ocean surface. The proportionality constant k, sometimes referred to as the piston velocity is generally tracer dependent. Note that the Dirichlet boundary condition is obtained asymptotically as k→∞. In other words, when the air-sea gas exchange process is sufficiently fast the Robin boundary condition is essentially the same as prescribing the concentration at the sea surface. A mixed boundary condition in which a linear combination of the flux and the surface tracer concentration is prescribed is often referred to as a Robin boundary condition (e.g. Haine, 2006). The interested reader is referred to Appendix D for the continuous formulation of the eigenvalue problem for the case of a Robin boundary condition applied at the sea surface. Here we present the discrete problem for the approach to equilibrium subject to a Robin boundary condition. Because air-sea gas exchange occurs over the entire sea surface we will only consider the case in which the Robin boundary condition is applied to the entire surface ocean and will not consider the localized-patch case as we did for the Dirichlet problem. The discrete tracer transport problem for a Robin boundary condition can be written in terms of a restoring timescale τ = z/k where z is the thickness of the top most layer of the ocean model. In matrix form we have where 1 is a vector of ones and where is a diagonal matrix with ones in the columns corresponding to surface grid points and zeros everywhere else. To obtain the discrete eigenmode problem for the Robin boundary condition we substitute the modal form c(t)=ve −λt in the homogeneous counterpart of Eq. (30) to obtain, where, as before, N=63 090, is the number of grid boxes in the ocean model.

Approach to equilibrium for the Robin problem
The solution of the Robin problem in Eq. (30) is mathematically similar to the Dirichlet problem considered in Sec. 3.1. If we enlarge our system to include the air-sea tracer exchange process, the Robin problem can be viewed of as a Dirichlet problem in which the tracer concentration is prescribed in the air overlying the sea surface instead of in the water at the sea surface. Not surprisingly, the solution to the Robin problem can be expressed as the sum of a forced plus homogeneous solution as was the case for the Dirichlet problem, where the leading 1 corresponds to the forced response due to the prescribed concentration of unity in the atmosphere, and where the remaining terms are the eigenmode expansion of the homogeneous solution. The eigenmodes in Eq. (33) that capture the approach to equilibrium are those of the enlarged system in Eq. (32) that includes the process of airsea exchange. The eigenmodes are therefore functions of the timescale τ with which the 50 m thick top layer of the model equilibrates with the atmosphere. As we have already mentioned the Dirichlet boundary condition is a limit case of the Robin boundary condition. We therefore expect that if the air-sea gas exchange time scale is small compared with the residence time of water in the mixed layer, the eigenmodes of the Robin problem will be similar to the eigenmodes of the Dirichlet problem. Figure 7 shows the e-folding decay timescale of the most slowly decaying eigenmode of the Robin problem as a function of the air-sea equilibration timescale for a 50 m thick mixed layer. As before the modes were computed using Matlab's eigs function. As τ →0, the e-folding timescale of the most slowly decaying eigenmode asymptotes to τ D =797 years, the efolding timescale of the most slowly decaying eigenmode of the Dirichlet problem with the surface patch being the whole sea surface. For values of τ shorter than 3 months, the efolding timescale of the most slowly decaying mode for the Robin problem is within 10% of the Dirichlet asymptote.
If on the other hand the air-sea equilibration timescale is sufficiently long, the reciprocal of τ is a small parameter, and the eigenvalue problem in Eq. (32) can be thought of as a perturbation to the eigenvalue problem for the Neumann problem discussed in Sect. 3.2. By having introduced the process of air-sea gas exchange, the tracer can escape to the atmosphere and is not conserved in the ocean as it was for the case of the Neumann problem. The process of air-sea tracer exchange therefore perturbs the null eigenmode of the Neumann problem into a decaying mode. A perturbation expansion in powers of 1/τ allows us to examine how the null eigenmode of the Neumann problem is modified by the process of air-sea gas exchange. Expanding the eigenvector and eigenvalue of the Robin probem in powers of 1/τ , substituting into Eq. 32 and equating like powers of 1/τ we recover to zeroth order the Neumann eigenvalue problem, Taking φ 0 and γ 0 to be the null mode of the Neumannproblem, we have φ 0 =1 and γ 0 =0. At first order we obtain, where we have used the fact that because of tracer conservation volume weighted sum of a column of the transport operator is zero, i.e. 1 WA = 0 . The numerator in Eq. (37) is equal to the volume for the top layer of the ocean model, and the denominator is equal to the total volume of the ocean. The large-τ asymptotic behavior of the most slowly decaying eigenmode is therefore where A is the surface area of the ocean and k is the airsea tracer-exchange rate. Equation (38) shows that when the air-sea exchange rate is sufficiently slow the tracer equilibration timescale for the ocean can be obtained simply by dividing the ocean volume by the rate of tracer input through the air-sea interface. The ocean, subject to a Robin boundary condition with a sufficiently slow air-sea tracer exchange rate, behaves essentially as a well-mixed box in terms of its equilibration timescale. The large τ asymptotic limit is indicated in Fig. 7 by the solid blue line. For values of τ greater than approximately 100 years the e-folding timescale of the most slowly decaying eigenmode is within 10% of the asymptotic value given in Eq. (38). For τ =100 years the efolding timescale of the most slowly decaying mode is 8225 years. For reference, the time for a 50 m thick mixed layer to come to equilibrium with the overlying air is on the order of a week for oxygen or chlorofluorocarbon and roughly 10 years for radiocarbon (Broecker and Peng, 1982). The long equilibration timescales obtained in the asymptotic limit of 1/τ →0 have of course more to do with the slow air-sea gas exchange process than with the circulation of the ocean. Figure 8 shows the spatial patterns associated with the most slowly decaying eigenmode for the Robin problem for several values of τ in the transition region between the small and large τ asymptotic limits. As τ increases, the spatial gradients of the most slowly decaying eigenmode gradually decrease, first in the deep ocean but eventually throughout the whole ocean. Eventually, for τ >500 years, the modal pattern is essentially constant. When the rate at which a tracer is introduced into the ocean is slow compared with the rate at which the ocean circulation can mix a tracer on global scales, the relative tracer disequilibrium between different parts of the ocean is necessarily small. Figure 9 shows the time evolution of a tracer subject to a Robin boundary condition with the tracer concentration in the atmosphere prescribed to unity for t>0. The air-sea equilibration timescale is set to τ =10 years, the value in the middle of the transition between the two asymptotic limits. The approach to equilibrium is slow, consistent with the correspondingly large e-folding timescale of 1690 years for the most slowly decaying eigenmode. The time for the tracer concentration to come withing 10% of its final equilibrium (t 90 ) is greater than 4000 years for both the deep Atlantic and Pacific oceans. Despite the long equilibration time, the lag-time for the Pacific to reach the same level of equilibration with the atmosphere as the Atlantic is only on the order of 1000 years. This is consistent with the eigenmode 24 F. Primeau and E. Deleersnijder: Tracer equilibrium pattern shown in Fig. 8. The Atlantic-to-Pacific gradients in the deep ocean are already substantially weakened for τ =10 years. The modal pattern for τ =10 years shows that the relative disequilibrium between the surface and deep ocean is expected to be much greater than the Atlantic-to-Pacific disequilibrium. Direct time-dependent simulation (not-shown) revealed that the time lag between the surface and deep ocean to reach 10% of their final equilibration value is on the order of 4000 years for τ =10 years, and only on the order of 2500 years for τ =1 month.
It is interesting to point out that the air-sea equilibration time for 14 C of τ ≈10 years lies right in the middle of the transition between Dirichlet and Neumann asymptotic limits in Fig. 7. This implies that ignoring the air-sea exchange process (small τ limit) or treating the ocean as a well mixed box (large τ limit) will result in significant underestimates in the equilibration time of 14 C. The importance of including both ocean circulation and air-sea exchange process for understanding the time evolution of 14 C is of course well known by the paleoceanographic community (e.g. .

Discussion
The type of surface boundary condition one should apply depends on the information available. The Dirichlet boundary condition is appropriate when the surface concentration is known. The Neumann boundary condition is appropriate when the flux is known. In the case where the tracer enters the ocean through air-sea gas exchange the information we have is often for the atmospheric concentration of the tracer. In that case a Robin boundary condition in which a linear combination of the flux and surface concentration is specified is appropriate.
Tracer simulations that use a prescribed tracer concentration as a surface boundary condition allow the tracer flux to be determined as part of the solution. This has the effect of making the flux dependent on the tracer concentration and on the ocean circulation. However, for a tracer such as δ 18 O, the flux of melt water into the ocean determines the initial surface concentration of δ 18 O and not the other way around.
The injection of a tracer by a process such as ice melt is consistent with prescribing a flux over a limited region of the surface ocean. Prescribing the concentration as is done by WH08, implies that there is some instantaneous feedback mechanism that keeps the surface δ 18 O concentration at some prescribed value as deeper waters with little or no tracer signature are mixed into the surface layer. It is difficult to imagine a mechanism that could achieve this. We expect no such feedback from the runoff of ice-melt from rivers and streams and while the atmosphere is a reservoir with a fast mixing timescale it is a rather small reservoir and is therefore not expected to have the ability to keep the δ 18 O content of surface waters constant during the period over which the ocean reaches equilibrium.
For a tracer such as 14 C that enters the ocean through air-sea gas exchange a Robin boundary condition is appropriate. To the extent that the air-sea gas exchange is sufficiently rapid, one might argue that a prescribed concentration boundary condition can be used to provide timescales useful for the interpretation of regional differences in tracer concentrations. However, for the specific case of 14 C the air-sea disequilibrium is sufficiently large that it needs to be taken into account  by either using a Robin-type boundary condition with an appropriate air-sea exchange rate or by prescribing a flux boundary condition in a well mixed atmospheric box coupled to the ocean.
While the Neumann boundary condition in which the flux of tracer is prescribed is the relevant one for understanding the equilibration time of a tracer such as δ 18 O, the Dirichlet boundary condition in which the concentration is specified is still extremely useful for understanding the age concept in marine modeling (e.g. Delhez et al., 1999;Holzer and Hall, 2000;Deleersnijder et al., 2001;Deleersnijder et al., 2002;Haine and Hall, 2002). As mentioned in the introduction, for the case of a stationary circulation the solution to the Dirichlet boundary condition can be interpreted as the cumulative distribution of times since the fluid elements in the ocean interior were last in contact with the surface patch. The time derivative of the cumulative distribution is often referred to as an age distribution and it gives important information about the distribution of times by which fluid elements are transported from the surface ocean to a point in the interior of the ocean via multiple pathways. However, in order for the age distribution to give specific information about the relative disequilibrium of a given tracer in the interior of the ocean it must be convolved with the time history of the surface concentration of the given tracer. Unless the concentration history of the tracer at the surface is known, the age distribution cannot be used to infer the relative disequilibrium of the tracer. This point is made explicit from the integral equation relating the Green function, G(r, t; ) for propagating a pulse of tracer through a surface patch at time t=0 and the Green function G(r, t; r , t ) for propagating a prescribed surface concentration at time t at point r inside , (The formulation of the governing equations that define G and G can be found in Holzer and Hall, 2000). Equation (39) shows how an initial pulse of flux through at time t=0 is propagated using G(r , t ; ) from back to a point r on at some time t ≥0 and then from r to a point r in the interior at time t using G(r, t; r , t ). The appearance of G(r , t ; ) on the right hand side of Eq. (39) shows that detailed information about the time evolution of the tracer concentration within is needed to obtain the interior tracer disequilibrium at time t. This information is not contained in G, the Green function for the Dirichlet problem. Holzer and Hall (2000) point out that the mean time since last contact with a patch , = ∞ 0 τ G(r, τ ; )dτ , tends to infinity as the size of shrinks to zero. In other words, as shrinks to zero it takes infinitely long on average for fluid elements to find their way from to r. It may therefore seem paradoxical that a tracer that is initially concentrated in a small patch can become well mixed in a finite time even as shrinks to zero. There is no paradox, however. It is the short τ part of the G distributions that carries the bulk of the tracer because as the area of the patch shrinks to zero, the initial concentration tends to infinity. In the case of the δ 18 O tracer for example, the concentration of δ 18 O depleted water is expected to have been extremely elevated at the mouth of the melt water rivers and streams relative to the well mixed asymptotic state. Therefore, when is small, it is the fast transport of waters with initially concentrated tracer values that carry the bulk of the δ 18 O signal and so bring the ocean to its well mixed state long before t 90 or is reached.

Conclusions
In this article we have used the eigenmodes of a simple three box model and of a three-dimensional ocean circulation model to characterise the timescale for a tracer to come to equilibrium. We have shown that the timescale to reach equilibrium is longer when one prescribes the concentration of tracer instead of the flux of tracer. The difference in equilibration time is especially pronounced when the patch over which the tracer is injected is small. For the Dirichlet problem the time to reach equilibrium tends to infinity as the area of the patch shrinks to zero. In contrast, the equilibration time of the Neumann problem becomes insensitive to the details of the injection as the patch size shrinks to zero. This result should not be surprising to anyone who has observed the dispersion of a dye in fluid experiment -the size of the syringe used to inject the dye quickly become irrelevant to the evolution of the dye. Prescribing the surface concentration to be a constant implies that the flux of tracer into the ocean depends on the concentration of this tracer in surface water. Such a boundary condition can be used for interpreting the time to uniformity of a tracer that enters the ocean through air-sea gas exchange, provided that the air-sea exchange rate is sufficiently rapid and that the Dirichlet boundary condition is applied to the entire ocean surface. Applying the Dirichlet boundary condition in patches of limited surface area breaks the superposition principle and makes it difficult to apply the resulting time scales for the interpretation of real tracers.
For the case of δ 18 O, the melting of ice sheets is causally independent of the concentration of δ 18 O in surface waters. The correct boundary condition for interpreting the relative disequilibrium between the deep Atlantic and Pacific oceans found in the δ 18 O record is the Neumann boundary condition for which the flux is prescribed. Our three-dimensional model results suggest that the relative disequilibrium between the modern day deep Atlantic and Pacific oceans is on the order of 1200 years or less. This tracer equilibrium timescale depends of course on the accuracy of our model's velocity field and eddy diffusivity coefficients. Given that the model resolution is very coarse one might wonder how accurately it can simulate the time-mean circulation. Early studies such as Cox, (1989) and England (1993) have shown that coarse-resolution models are capable of representing the global-scale water masses suggesting that at least on the largest scales the model circulation is reasonably accurate. Given that the most slowly decaying eigenmodes are associated with the largest spacial scales we do not expect that our result should be very sensitive to the model resolution. Of course, a realistic simulation of the dispersion of the δ 18 O signal from the melting of the ice sheets needs to take into account the changes in ocean currents that are expected to occur as a result of the change in density due to the fresh-water pulse. Without taking these effects into account we can still question the "simpler explanation" offered by WH08 that a steady modern-day circulation is sufficient to account for the apparent 3900 year lag between the times for the deep Pacific and Atlantic oceans to reach equilibrium.

Formulation of the governing equations for the 3-box model
The volumes of the boxes are denoted by V 1 , V 2 and V 3 , and the corresponding tracer concentrations are denoted by c≡[c 1 c 2 c 3 ] . As indicated in Fig. 2, the index "3" is associated with the deep ocean box, while the indices "1" and "2" are associated with the two surface ocean boxes. The geometry of the boxes, as shown in Fig. 2a, is such that S denotes the area of the interface separating the deep box from the surface boxes as well as the area of the ocean-atmosphere interface. The limited area of the interface separating boxes 1 and 3 is then αS with 0≤α≤1, and the area of the interface between boxes 2 and 3 is (1−α)S. To keep the algebraic complexity to a minimum it is useful to consider the case where the thickness of the surface boxes is much less than that of the deep box. We therefore introduce the following small parameter The net flux of tracer from box i to box j is denoted by φ i,j as indicated in Fig. 2b, and the net flux of tracer from the atmosphere into surface box "1" is denoted by φ a,1 . The flux of tracer from the atmosphere to box "2" is assumed to be zero. For simplicity and since is assumed small we will neglect the tracer exchange between the two surface boxes, more precisely, we set φ 1,2 =0. We parameterize the fluxes between the ocean boxes as the product of (i) the area between the boxes, (ii) a constant velocity scale v and (iii) the difference in tracer concentration between the boxes: The tracer budget for each of the three boxes yields the following governing equations Combining the flux relations in Eq. (A2) with the differential equations in Eq. (A3) we obtain where we have introduced the time-scale It is convenient to rewrite this system of equations in matrix form dc dt + Ac = s, Appendix B

Continuous problem for the case of a Dirichlet boundary condition
The mathematical formulation of the tracer initial value problem with concentration prescribed over a fixed patch (i.e. the problem considered by WH08) is given by ∂ ∂t C(t, r) = −∇ · [u − K · ∇] C(t, r), subject to the boundary conditions n · (K · ∇C(t, r)) = 0 for r on the boundary ∪ 1 \ 0 , C(t, r) = 0, for t ≤ 0 1, for t > 0 on r ∈ 0 , where n : is a unit vector normal to the basin boundary, : is the bottom and sides of the ocean basin, 0 : is a surface patch at the air-sea interface, 1 : is the total air-sea interface, 1 \ 0 : is the rest of the air-sea interface excluding the 0 patch.
The eigenfunctions, ψ n (r), of the Dirichlet problem are the solutions to the following eigenvalue problem      (−1/τ n + iω n )ψ n (r) = −∇ · [u − K · ∇]ψ n (r), ψ n (r) = 0 for r on the boundary 0 , n · (K · ∇ψ n (r)) = 0 for r on the boundary 1 \ 0 , where the eigenmodes are ordered such that their decay rate increases with increasing n. Note that because of the absorbing boundary condition over 0 and the absence of any tracer source terms, the eigenmodes, e λ n t ψ n (r), n=1, 2, · · · , ∞, must be decaying functions of time, implying that τ n >0 for all n. : is the bottom and sides of the basin boundary, 1 : is the total air-sea interface, 0 : is a surface patch on the air-sea interface, δ(t) : is the Dirac-delta function, V tot : is the total volume of the ocean, A 0 : is the area of the 0 patch.
The eigenfunctions, ψ n (r), of the Neumann problem are the solutions of the following eigenvalue problem (−1/τ n + iω n )ψ n (r) = −∇ · [u − K · ∇]ψ n (r), n · (K · ∇ψ n (r)) = 0 for r on the boundary ∪ 1 , where the eigenvalues, λ n =(−1/τ n +iω n ), n=0, 1, · · · , ∞, are ordered such that the e-folding decay timescale, τ n , decreases with increasing n. As already mentioned in Sect. 2.1, the zero eigenvalue is a manifestation of the conservation of tracer. That the eigenfunction corresponding to the zero eigenvalue is constant is consistent with the physical fact that in the absence of source or sinks, a uniform tracer distribution must be in steady state. Also, because diffusion acts to decrease the difference between the domain maximum and minimum tracer concentrations, all the eigenmodes except for the one with the zero eigenvalue must be decaying modes, (i.e. τ n >0). Conservation of tracer mass then requires that all eigenfunctions with n≥1 integrate to zero. This last fact implies that any initial tracer distribution that does not integrate to zero must necessarily project onto ψ 0 .