The Ocean’s Memory of the Atmosphere: Residence-Time and Ventilation-Rate Distributions of Water Masses

A conceptually new approach to diagnosing tracer-independent ventilation rates is developed. Tracer Green functions are exploited to partition ventilation rates according to the ventilated fluid’s residence time in the ocean interior and according to where this fluid enters and exits the interior. In the presence of mixing by mesoscale eddies, which are reasonably represented by diffusion, ventilation rates for overlapping entry and exit regions cannot meaningfully be characterized by a single rate. It is a physical consequence of diffusive transport that fluid elements that spend an infinitesimally short time in the interior cause singularly large ventilation rates for overlapping entry and exit regions. Therefore, ventilation must generally be characterized by a ventilation-rate distribution, (cid:1) , partitioned according to the time that the ventilated fluid spends in the interior between successive surface contacts. An offline forward and adjoint time-averaged OGCM is used to illustrate the rich detail that (cid:1) and the closely related probability density function of residence times R provide on the way the ocean communicates with the surface. These diagnostics quantify the relative importance of various surface regions for ventilating the interior ocean by either exposing old water masses to the atmosphere or by forming newly ventilated ones. The model results suggest that the Southern Ocean plays a dominant role in ventilating the ocean, both as a region where new waters are ventilated into the interior and where old waters are first reexposed to the atmosphere.


Introduction
The rate of water-mass ventilation and its spatial distribution are fundamental to understanding the role of the ocean in climate. Ventilation fluxes are immediately relevant to 1) determining the uptake of anthropogenic carbon dioxide by the ocean (e.g., Wallace 2001;Gruber et al. 1996;Hall et al. 2004), 2) constraining high-latitude deep-water formation rates and possibly thermohaline variability using transient and radioactive tracers (e.g., Broecker et al. 1999;Orsi et al. 2001), and 3) studying the ventilation and circulation of the subtropical thermocline (e.g., Jenkins 1998;Robbins et al. 2000).
The surface fluxes of heat and moisture that drive water-mass transformations are difficult to measure. Consequently, formation rate estimates are very uncertain (Large and Nurser 2001). This is particularly true for deep-water masses for which the measurement difficulty is compounded by entrainment: "The very definition of bottom-water formation rates is vague owing to the importance of entrainment in forming water mass properties" (Munk and Wunsch 1998). The same point is emphasized by Wunsch and Ferrari (2004). Even less is known about the pathways and mechanisms by which fluid returns to the surface.
Tracers offer an alternative method for determining ventilation rates because transient and radioactive species such as tritium, chlorofluorocarbons (CFCs), and 14 C label fluid elements with their near-surface concen-trations. However, deriving water-mass formation rates from tracer measurements is complicated by uncertain estimates of eddy diffusion. Diffusion plays a crucial role in transport by mixing waters that last equilibrated with the surface at different times and places. While it is common to treat transport as purely advective when computing ventilation rates, the neglect of eddy diffusion leads to substantial errors (Hall et al. 2004). Moreover, we will argue here that diffusion prevents quantification of a ventilation rate in terms of just a single number.
To interpret tracer measurements in terms of watermass formation rates, one must use models of the ventilation process that make assumptions (often implicit) about the importance of diffusive fluxes. Because the assumptions made in different studies are not necessarily consistent with each other or with the real ocean, estimates of water-mass formation based on different tracers can yield different results. For example, using CFC and 14 C data, very different estimates of deepwater formation rates have been obtained with box models by Broecker et al. (1999) versus Orsi et al. (2001).
Because box models cannot fully capture turbulent mixing, entrainment, and eddy diffusion, it is not clear to what extent a ventilation rate derived from a box model using a particular tracer can represent a tracerindependent ventilation rate that is a property of only the fluid flow. Properly defined ventilation rates should measure the rate at which water recently exposed to the atmosphere enters the ocean interior and the rate at which interior water becomes reexposed to the atmosphere, independent of any tracers used to determine these rates. Once the ventilation rates are known, they can be used to calculate the uptake and reexposure of any trace species of interest.
In this paper, we develop tracer-independent diagnostics of the rate of ventilation into and out of the ocean interior. In the process, we introduce a new way of thinking about ocean ventilation that recognizes from the outset that diffusive mixing will give fluid elements fluxed into the ocean a broad distribution of interior residence times. We begin by considering ventilation of the entire ocean domain through the global ocean surface. The flux of water ⌽() that has a residence time in the interior of or longer is shown to be equivalent to the global inventory of water binned according to age (i.e., time since last surface contact). In the continuum limit of advective-diffusive flow, the ventilation rate for fluid regardless of its residence time, ⌽(0), becomes infinite; that is, the ventilation rate is overwhelmed by fluid that spends vanishingly small time in the interior before reventilating. Therefore, a single ventilation rate generally does not suffice to quantify any water-mass formation or reexposure rates. The entire distribution ⌽() must be considered. The same issue also arises for estimating the one-way (or "gross") flux of stratosphere-troposphere exchange as shown by Hall and Holzer (2003). While there are precise parallels to the work of Hall and Holzer (2003) for a globally integrated ocean domain (see also Hall et al. 2005, manuscript submitted to J. Phys. Oceanogr.), our work here substantially further develops and generalizes the basic ideas.
We utilize a description of ocean transport in terms of boundary-propagator Green functions (Holzer and Hall 2000;Haine and Hall 2002;Primeau 2005), to define a ventilation-rate density distribution, (, ⍀ i , ⍀ f ), that partitions the flux with which new water masses are formed or older water masses are reexposed to the atmosphere, according to the water-mass residence time and places of formation (⍀ i ) and reexposure (⍀ f ). We show that, if formation and reexposure regions overlap, then (, ⍀ i , ⍀ f ) has a singularity at ϭ 0. This singularity is a physical consequence of diffusion and not an artifact of our formulation. The singularity implies that, if water is considered for which entry and exit regions overlap, then ventilation cannot be quantified by a single number unless some additional conditioning is imposed, such as a minimum time that the water must spend in the interior. Conditioning ventilation rates on nonoverlapping entry and exit regions removes the singularity because such regions cannot be connected by water with infinitesimal residence time. We also show that is very closely related to the probability density function (pdf) of residence times, R . Unlike , the pdf R has a physically meaningful normalization even if entry and exit regions overlap.
Throughout, we use the steady-state offline ocean model of Primeau (2005) to illustrate what our diagnostic framework reveals about the ocean circulation. To the extent that the model circulation is in agreement with the real circulation, our analysis paints a detailed and quantitative picture of what fraction of the ocean is involved in connecting different surface regions and how this fraction is further partitioned according to residence time. The structure of the residence time pdf's points to the underlying transport pathways. Overall, we give a novel and quantitative characterization of the ocean's surface-to-surface transport that is broadly consistent with our current knowledge of the large-scale circulation. A brief discussion of some published estimates of water-mass formation rates obtained from transient tracer measurements is presented in light of our theory.

Residence-time-partitioned ventilation rate: The basic idea
Depending on the ocean currents and eddy mixing processes that act on a parcel of fluid after it has been ventilated at the surface, some of the parcel's fluid elements are reexposed to the atmosphere after only a very short residence time below the surface, while others are transported far from the formation site before first being reventilated at the surface. The residence time of fluid elements, defined here as the time between successive surface contacts, is a useful transport diagnostic for analyzing the air-sea exchange of transient atmospheric gases because the time between successive ventilations determines the extent of their disequilibrium with the atmosphere.
The ocean is a sink for atmospheric constituents where they are ventilated into the ocean and a source of these constituents where they are ventilated out of the ocean. Therefore, it is natural to ask at what rate a given water mass that last ventilated into the ocean through some surface patch ⍀ i , and that has resided in the ocean interior in a given residence-time interval, will be reventilated or exposed to the atmosphere through another patch ⍀ f . Symmetrically, one would also like to know the rate at which the water mass that is ventilating out of the ocean at ⍀ f formed by ventilating into the ocean at ⍀ i .

a. Identifying and tracking ventilated fluid
To label fluid elements based on where and when they ventilate at the sea surface we use the boundary propagator, G (r, t | ⍀ i , t i ), or its counterpart for the timereversed adjoint flow, G (r, t| ⍀ f , t f ) (Holzer and Hall 2000;Primeau 2005). 1 The forward boundary propagator G labels fluid elements as they are ventilated into the ocean at surface patch ⍀ i and time t i . The equation of motion for G is subject to the surface boundary condition and where is the mass flux of G for fluid velocity u and eddy diffusion tensor K. At boundaries other than the surface, J · n ϭ 0, where n is a unit vector normal to the boundary.
The propagator G (r, t | ⍀ f , t f ) for the time-reversed adjoint flow labels fluid elements as they become exposed to the atmosphere at time t f at surface patch ⍀ f . The equation of motion for G is the same as (1), but with the delta function ␦(t Ϫ t f ) applied on ⍀ f and with where ١ · u ϭ 0 is assumed and K T denotes the transpose of K. The equation for G is integrated backward in time to where fluid elements labeled by G were last ventilated.
It is sometimes convenient to refer to the collection of fluid elements labeled by G as a "newly ventilated water mass" carrying the surface properties at patch ⍀ i and time t i . Similarly, one may call the collection of fluid elements labeled by G an "old water mass" that becomes reexposed to the atmosphere at patch ⍀ f and time t f .

b. Global inventory and ventilation
It is helpful to introduce the basic idea behind a residence-time-partitioned ventilation rate in its simplest form. To this end, this section specializes to the case of a stationary flow and ventilation through exposure to the atmosphere over the entire global sea surface. There is no need to distinguish between ventilation into or out of the entire ocean surface since both rates are equal in steady state. Local ventilation rates for nonsteady flow are developed in the next section.
For steady flow, G (r, t i ϩ | ⍀ i , t i ) depends on time only through the transit time, or "age," ϭ t Ϫ t i . Hence, for steady flow and ⍀ i ϭ ⍀, the entire ocean surface, we define G (r, ) ϵ G (r, t i ϩ | ⍀, t i ). In this case, G is the transit time distribution of times since last contact with the surface anywhere. The residence time of a fluid element "born" on ⍀ at time t i when it had last contact with the surface is then given by its age at the time of subsequent surface contact when it loses its surface label by virtue of the boundary conditions on G and "dies." Because G (r, ) d can be interpreted as the mass fraction at r of water with an age in the interval (, ϩ d), the mass flux of water labeled by G onto ⍀ is the rate with which water with a residence time in the interval (, ϩ d) is being reexposed to the atmosphere. In steady state this is also the rate at which newly ventilated waters are formed. 1 In the notation of Primeau (2005) is the joint pdf at r for last-passage time |p from ⍀ i , and ͐ ⍀ f dr f F(r f , fp ϭ t Ϫ t f | r, t)ϵG(r, t | ⍀ f , t f ) is the joint pdf at r for first-passage time fp ⍀ f . In the notation of Holzer and Hall (2000) and Hall and Holzer (2003), G † ϵ G, with the dagger instead of the tilde labeling solutions for the time-reversed adjoint flow. We changed the notation to avoid potentially confusing G † with the adjoint of G.
The ventilation rate for water that has resided in the interior for or longer is obtained by integrating the flux of G normal to ⍀ starting at time t ϭ after the water mass labeled by G was formed at t ϭ 0: where da is the ocean-to-atmosphere directed normal surface area element. Except for the sink by ventilation at the surface, G (r, ) evolves as a conservative tracer for Ͼ 0. Therefore, all the G tracer in the ocean at time t ϭ will ultimately flux out through ⍀ from time onward. Consequently, ⌽(), which is the accumulated flux of G leaving the ocean from onward, must equal the total tracer mass in the ocean at time : where V is the ocean volume. Thus, ⌽() has two complementary interpretations: (i) ⌽() is the cumulative ventilation rate of water that has resided in the ocean interior for or longer and, because (5) shows ⌽() to be the domain-integrated distribution of transit times ("ages") since last surface contact, we have that (ii) ⌽()d is the mass of the ocean whose age is in the interval (, ϩ d). We can, therefore, also think of ⌽() as the global inventory of water binned according to its age (see Fig. 1). A ventilation rate partitioned according to the residence time of the fluid being fluxed through the surface is given by the derivative of the cumulative flux: Here () characterizes the ocean's ventilation rate as a density distribution such that ()d gives the mass flux of newly ventilated water that will reside in the ocean interior for a time ∈ (, ϩ d) before being reventilated.
The equivalence between (4) and (5) may formally be obtained by integrating the equation of motion of G and applying Gauss's law, which was also derived in the context of stratosphere-troposphere exchange by Hall and Holzer (2003). However, the construction of a ventilation-rate density distribution, (), and the development that follows is new and has no existing precedents in the atmospheric or oceanographic literature to the best of our knowledge.

c. Connection with the residence-time pdf
The ventilation-rate density distribution () has a natural connection with R ()d, the instantaneous mass fraction of the ocean that has a residence time ∈ (, ϩ d). Figure 1 shows ⌽() schematically as a histogram of the oceanic population of fluid elements, or simply particles, binned according to age. As time advances, particles move to the right from one age bin to the next as they get older, unless they die through contact with the surface. The age of a particle at death is its residence time. The mass difference, ⌬M() ϵ Ϫ⌬⌬⌽(), by which an age bin at is shorter than its neighbor to the left, is the mass of particles that die with a residence time ∈ (, ϩ ⌬). For a steady ocean, every bin from 0 to contains a mass ⌬M that will ventilate/die when it reaches the bin centered on . Therefore, the total mass of water in the ocean at any instant that will have residence time (age at death) is given by n⌬M, where n ϭ /⌬ is the number of bins between age 0 and . Because by definition, R ()⌬ ϭ n⌬M/M, where M is the total ocean mass, we have R ⌬ ϭ Ϫ⌬⌽/M or, equivalently, taking the continuum limit ⌬ → 0 and using Ϫd⌽/d ϭ , we have Equation (7) simply states that, if the steady flux of fluid whose residence time lies in the interval (, ϩ d) acts for a time , the mass of the ocean present at any given time whose residence time falls in the same interval will be flushed out. Written as ϭ MR ()/(), this is a "spectral" version (i.e., a version holding sepa-FIG. 1. Schematic of the cumulative ventilation-rate distribution, ⌽(), which shows ⌽() as a histogram that bins the mass of the ocean in each age class of width ⌬. The dark shading indicates the mass of water in the age bin ( Ϫ ⌬, ) that will be exposed to the atmosphere during the next time interval, (, ϩ ⌬). Thus, this water has a residence time in the interval (, ϩ ⌬). At any given time the mass fraction of the ocean that will have a residence time in the interval (, ϩ ⌬) is given by R ()⌬. This fraction is indicated by the shaded rectangle of area Ϫ⌬⌽/M, where M is the total mass of the ocean. rately for each residence-time interval, or "class") of the familiar steady-state result that the residence time is given by the total mass divided by the rate of input.

d. Singularities of ⌽, R , and
The domain-integrated boundary propagator, ͐ V d 3 r(r)G (r, ), which we showed to equal ⌽(), has a known Ϫ1/2 singularity at ϭ 0 (Holzer and Hall 2000;Hall and Holzer 2003). Correspondingly, as → 0, has a Ϫ3/2 singularity [cf. (6)] and R a Ϫ1/2 singularity [cf. (7)]. We emphasize that these singularities are not artifacts of the way our diagnostics are defined, but a physical consequence of diffusion. If diffusing fluid elements are regarded on small scales as random walkers, then the singularity comes from the fact that such walkers make an ever-increasing number of contacts with the surface during any finite time interval as the continuum limit is approached at constant diffusivity (Hall and Holzer 2003).
Note that the Ϫ3/2 divergence of () is not integrable, which has important implications for estimating ventilation rates from tracer data (e.g., Broecker et al. 1998;Broecker et al. 1999;Orsi et al. 2001). Explicitly, as 0 → 0. Thus, the rate at which a newly ventilated water mass is formed, regardless of residence time, that is, 0 ϭ 0, cannot meaningfully be estimated from tracer data because it is infinite. While ⌽( 0 ) is finite for 0 Ͼ 0, there is generally no natural choice for the cutoff, 0 , and the divergence of ⌽ implies high sensitivity on the value of 0 . Because R () is a proper pdf of residence time , with finite normalization (and physically meaningful moments; see section 4a), the Ϫ1/2 divergence of R is integrable. The physical difference between the distributions R and is that they partition different particle populations. The distribution R () bins the finite mass of the ocean according to the birth-to-death lifetime of the composing water particles at an instant in time, ensuring a properly normalized pdf. By contrast, the residence-time-partitioned ventilation rate (), acting for a clock-time interval ⌬t, that is, ()⌬t, bins according to the mass of fluid that makes surface contact during ⌬t. The crucial difference between R and ⌬t is that individual fluid elements with Ͻ ⌬t contribute just once to R , but more than once to ⌬t. Particles with Ͻ ⌬t make repeated surface contact during ⌬t, whence they die and are instantly reborn under a new identity. Indeed, in the continuum limit with diffusion, particles with approaching zero make an everincreasing number of contacts with the surface during any finite ⌬t so that the size of the population partitioned by ()⌬t with ∈ ( 0 , ϱ) grows unbounded as 0 → 0.

e. Illustration with global ocean model
As we develop our framework, it is useful to illustrate what the diagnostics reveal about the ocean. To this end we use the off-line ocean transport model developed by Primeau (2005). This model is based on the time-averaged velocity and eddy-diffusivity tensor fields from a dynamical simulation with a full OGCM. The OGCM used is a version of the ocean component of the climate model of the Canadian Centre for Climate Modeling and Analysis, itself based on version 1.3 of the National Center for Atmospheric Research's (NCAR) Climate System Model (CSM) Ocean Model (Pacanowski et al. 1993). It uses the k-profile paramterization (KPP; Large et al. 1994) vertical mixing scheme and the GM (Gent and McWilliams 1990) isopycnal eddy-mixing scheme. The off-line model (like its parent OGCM) uses second-order centered differences on an ϳ3.75°ϫ 3.75°grid with 29 levels ranging in thickness from 50 m near the surface to 300 m for the deepest level.
Time-evolving forward and adjoint boundary propagators were computed by exponentiating the steadystate transport operator over finite time steps using the Expokit software of Sidje (1998). Moments of pdfs were computed by direct matrix inversion of their steady-state equations using the multifrontal approach of the multifrontal massively parallel sparse (MUMPS) direct solver (Amestoy et al. 2001).
Log-log plots of the model-derived global residencetime pdf R () and the cumulative ventilation-rate distribution ⌽() are shown in Fig. 2. The corresponding ventilation-rate density distribution () is shown in Fig. 3. Here ⌽(), which is also the global age inventory, shows remarkably little structure with a near-powerlaw scaling from ϳ10 to ϳ600 yr followed by a rapid falloff. From a globally integrated perspective, the time-averaged model circulation appears to behave very similarly to a simple one-dimensional diffusive reservoir. The global residence-time pdf R () and ventilation-rate density distribution () ϭ Ϫd⌽()/d show more structure, hinting at the importance of advective transport that, as we will see, is only fully revealed when we further partition according to where water is ventilated into, and out of, the ocean interior.
With decreasing , we see that R (), ⌽(), and () tend to diverge. The expected small-power laws are indicated in Figs. 2 and 3. Deviations from these scalings for Շ 10 yr can be attributed to the effects of the model's finite spatial resolution. But even at finite spatial resolution, the physical continuum-limit divergence manifests itself sufficiently strongly to make it clear that attempting to integrate over all is not physically meaningful.

Ventilation-rate distributions for time-dependent flow and specified formation and reexposure regions
We now relax the assumption of stationary flow and further partition the ventilation rates according to specified surface patches ⍀ i and ⍀ f where fluid elements are ventilated into and out of the ocean.

a. Generalization of the ventilation-rate density distribution
In a time-dependent flow, the mass fraction at time t that was formed on ⍀ i a time ago , with which waters are reexposed to the atmosphere Conservation of mass guarantees that the outward flux over ⍀ f at t f of water formed on ⍀ i at time t i ϭ t f Ϫ must be equal to the flux of water into the ocean on ⍀ i at time t i that is destined to be reventilated on ⍀ f at time t f . Thus, the residence-time-partitioned ventilation-rate distribution with which water forms on ⍀ i at t i that will have a residence time in the interval (, ϩ d) before reventilating on ⍀ f is given by the symmetry where t f Ϫ t i ϭ . In our notation, a down arrow (up arrow) denotes a surface flux into (out of) the interior. The arguments on the left of the semicolon refer to the time and place at which the flux is occurring; the arguments on the right of the semicolon refer to the identity of the fluid being fluxed, that is, its residence time and either its place of origin or place of first reexposure to the atmosphere. It is sometimes of conceptual and computational advantage to use the propagator of the adjoint equation, G, to obtain the ventilation fluxes. While in (9) one computes the instantaneous mass flux onto ⍀ f at t f of fluid elements labeled at (⍀ i , t i ), we can also label fluid elements at (⍀ f , t f ) and trace them through the timereversed adjoint flow to obtain The entire-surface-to-entire-surface ventilation-rate density distribution (). To guide the eye, the thin dashed line indicates a Ϫ3/2 power law.
FIG. 2. The entire-surface-to-entire-surface residence-time pdf R () (thick solid line) and the corresponding cumulative ventilation-rate distribution ⌽() (thick dashed line), which is the mass flux of water that will remain in the interior for at least time . Equivalently, ⌽() is the global inventory of water binned according to its age (transit time since last surface contact). To guide the eye, the thin dashed line indicates a Ϫ1/2 power law. Both R () and ⌽() are seen to diverge with decreasing over the range plotted. Neither R nor ⌽ obey a true power law. (If they did, it would have to be identical for both because R and ⌽ are related through R ϭ Ϫ‫ץ‬⌽/‫).ץ‬ Note that an entry patch for the forward flow becomes an exit patch for the time-reversed adjoint flow.
One can now straightforwardly derive the timedependent and patch-to-patch generalized versions of the cumulative flux ⌽. When either ⍀ f (or ⍀ i ) are the entire ocean surface, one can again show that ⌽ corresponds to a domain integral of G (or G ). Likewise, ↓ and ↑ are closely related to the residence-time pdf R . The patch-to-patch version, R (, ⍀ i , ⍀ f ; t), represents the joint pdf of residence time, , and entry and exit patches, ⍀ i and ⍀ f , so that R (, ⍀ i , ⍀ f ; t)d represents the mass fraction of the entire ocean at time t that "entered" through ⍀ i and will exit through ⍀ f after residing in the ocean interior for a time ∈ (, ϩ d).
The relations defining ↑ , ↓ , ⌽ ↑ , ⌽ ↓ , and R, in terms of both J and J , are collected in Table 1. The small-singularities of the distributions are derived in appendix A. If the flow is stationary, ↓ and ↑ are time independent so that we can define (, . Stationarity allows us to make the patch-to-patch generalized statement of (7), which for each residence-time class the residence time is the total mass divided by the rate of input:

b. Model-derived maps of ventilation
The residence-time-partitioned, local formation, and reexposure ventilation fluxes are given by the integrands of the surface integrals (T3) and (T2), respectively. The local ventilation fluxes regardless of residence time can then be obtained by integrating over all residence times. Thus, Fig. 5a (described below) shows ͐ ϱ 0 d J (r, 0| ⍀ f , ) · n, that is, maps of the residence-time-integrated formation flux of newly ventilated water spatially partitioned according to the surface patch, ⍀ f , where it will first be reexposed to the atmosphere. Figure 5b (described below) shows ͐ ϱ 0 d J(r, 0| ⍀ i , Ϫ ) · n, that is, maps of the residence-timeintegrated exposure flux of old water partitioned according to the surface patch, ⍀ i , where the water was formed. Stationarity allows us to choose t i ϭ 0 for J , and separately t f ϭ 0 for J. Our geographical choice for partitioning the ocean surface into patches is shown in Fig. 4.
The ventilation flux patterns are consistent with our understanding of the large-scale ocean circulation. For example, high exposure rates are observed in the upwelling regions of the Southern Ocean and eastern equatorial Pacific, and low exposure rates are observed in the subtropical gyres where the Ekman transport is convergent. Waters formed in the high-latitude North Atlantic (⍀ 1 ) and the Southern Ocean (⍀ 5 ), regions that dominate deep-water formation, contribute to high exposure fluxes in most upwelling regions, as expected. Other formation regions are also important contributors to upwelling in specific regions. For example, waters last exposed to the atmosphere in the TABLE 1. Table of quantities related to the formation of newly ventilated water masses (right column) and to the exposure to the atmosphere of old water masses (left column). The ventilation-rate density distribution ↓ is the mass flux through ⍀ i at time t i of newly ventilated water partitioned according to that water's residence time, , and location of ultimate reexposure to the atmosphere, ⍀ f . The ventilation-rate density distribution ↑ is the mass flux onto ⍀ f at time t f of old water partitioned according to that water's residence time and formation location, ⍀ i . The cumulative ventilation-rate distributions ⌽ ↓ and ⌽ ↑ are the corresponding mass fluxes with a residence time, , or longer. If ⍀ f ϭ ⍀ (the entire ocean surface), then ⌽ ↓ can also be thought of as the volume-integrated propagator, G, of surface boundary conditions (T13) and, if ⍀ i ϭ ⍀, ⌽ ↑ is the volume-integrated boundary propagator, G, for the time-reversed adjoint flow (T14). Here, R is the residence-time pdf, which is the mass fraction of the ocean at time t with residence time in the interval (, ϩ d) and entry and exit regions ⍀ i and ⍀ f . (In addition, M is the mass of the ocean, is the fluid density, and J and J are the fluxes of G and G.)

Quantities related to ventilation:
Into the ocean Out of the ocean Comparing Figs. 5a and 5b, we see that the subtropical gyres, which are evident in Fig. 5a as regions of elevated formation rate, are in Fig. 5b regions of low exposure rate. Conversely, the equatorial regions that appear as regions of elevated exposure rate in Fig. 5b are revealed as regions of low formation rate in Fig. 5a. This correspondence is a manifestation of advective fluxes in those regions. Such a correspondence does not exist for the Southern Ocean, which shows elevated rates of both formation and exposure. This is consistent with the importance of eddy diffusive fluxes associated with mixing and turbulent entrainment for Southern Ocean water-mass formation and exposure.
Note that in Figs. 5a and 5b we have avoided contouring the residence-time-integrated ventilation fluxes where the entry and exit regions overlap. The total flux for overlapping entry and exit patches is singular because of the overwhelming contribution from fluid el-ements with vanishingly small residence times in the interior (appendix A). The singular contribution to the total flux from fluid elements with → 0 masks the potentially significant contributions from fluid elements with finite residence times. For example, calculations based on the residence-time pdf R presented in section 4b show that 27.8% of the total ocean volume is ventilated by elements that enter and exit the ocean interior through the Southern Ocean (⍀ 5 ).
Patch and residence-time-integrated ventilation rates are organized in the table in Fig. 6. This table is simultaneously a table of formation and exposure rates because the ⍀ i → ⍀ f formation rate on ⍀ i is equal to the ⍀ i → ⍀ f exposure rate on ⍀ f . As discussed in appendix A, if the entry and exit patches do not overlap or share a border, patch-to-patch ventilation rates are well defined and physically meaningful (boldface numbers in Fig. 6). We see, for example, that ϳ6.5 Sv (Sv ϵ 10 6 m Ϫ3 s Ϫ1 ) of newly ventilated water formed in the highlatitude North Atlantic (⍀ 1 ) will be reexposed to the atmosphere in the Southern Ocean (⍀ 5 ).
Ventilation rates for patch pairs that are not physically meaningful (smaller font) are finite only because of the spatial discretization of the continuous equa-  (Gurney et al. 2000) and were also utilized in a study of ocean transport by Primeau (2005). tions. We show them here to emphasize that, even in a coarse-resolution model, ventilation rates in excess of 10 4 Sv can be obtained for coincident entry and exit regions (diagonal entries).

Steady-state residence-time pdfs and ventilation volumes a. Moments of R and ventilation volumes
Because R is a proper pdf, we can summarize some of its key features by calculating its moments: The temporal norm, r 0 (⍀ i , ⍀ f ), represents the mass fraction of the ocean at any given time that has formed on ⍀ i and will be reexposed on ⍀ f . The mean residence time, (⍀ i , ⍀ f ), of fluid formed on ⍀ i and reexposed on ⍀ f is given by (⍀ i , ⍀ f ) ϭ r 1 (⍀ i , ⍀ f )/r 0 (⍀ i , ⍀ f ), which provides a convenient time scale for advective-diffusive ⍀ i → ⍀ f transport. The ⍀ i → ⍀ f pathways with short connect parts of the oceans that are ventilated rapidly, whereas pathways with long connect parts of the ocean that ventilate slowly. The variance of , given by (r 2 Ϫ r 2 1 )/r 2 0 , is a bulk measure of the dispersion that occurs en route from ⍀ i to ⍀ f . For stationary flow, it turns out to be of great computational advantage to express the moments of R as fluxes of the moments of G. It follows immediately from (3), (T1), and (T5) that where J ␥m is the flux of the mth moment of G and J ␥m is the mth moment of G. In particular, when ⍀ i is the FIG. 6. Residence-time-and patch-integrated water-mass exposure/formation rates, ͐ ϱ 0 d(⍀ f , | ⍀ i , 0) (Sv), for all combinations of entry and exit patches. Rates below 0.01 Sv have been omitted for clarity. The numbers in small type would be infinite in the continuum limit and are therefore ill-defined. We show them here merely to emphasize that a model with finite resolution already gives huge diagonal terms (coinciding entry and exit patches). The expected singularity when ⍀ i and ⍀ f share an edge is only logarithmic and is not apparent at the resolution of the model (corresponding rates below 10 Sv have been omitted for clarity). entire ocean surface, ⍀, the first moment of G is the mean time since last ⍀ contact, also known as the ideal age (Holzer and Hall 2000;Haine and Hall 2002;Khatiwala et al. 2001), so that r 0 (⍀, ⍀ f ) is the flux through ⍀ f of the ideal age tracer. Similarly, r 0 (⍀ i , ⍀) is the flux through ⍀ i of the mean time to first ⍀ contact. [These particular relationships were used, but not derived, in (17) and (18) of Primeau (2005).] Relationship (15) allows us to compute spatially detailed properties of R with only a few tracers. The advantages are compounded if one can-as we have done here-recursively compute the moments of G and G by directly inverting the transport operator and its adjoint. For our coarse-resolution global ocean model, computational times are reduced by several orders of magnitude. If ⍀ f or ⍀ i is the entire ocean surface, ⍀, the moments of R are particularly simply related to the temporal moments of the domain-integrated G and or G : is twice the domainaveraged mean time since last ⍀ i contact, that is, twice ⌫(⍀ i ) ϵ ⌫ 1 (⍀ i )/ ⌫ 0 (⍀ i ). In the case in which ⍀ i ϭ ⍀, ⌫(⍀) is the domain average of the ideal age. Similarly, the mean ⍀ → ⍀ f residence time is twice ⌫ (⍀ f ) ϵ ⌫ 1 (⍀ f )/⌫ 0 (⍀ f ), the mean time to first ⍀ f contact. FIG. 7. Plots of the residence-time pdf, R (, ⍀ i , ⍀ f ) (yr Ϫ1 ), vs residence time (yr). The scale for is linear from 0 to 1000 yr with ticks every 200 yr, while the scale for R is linear from 0 to 5 ϫ 10 Ϫ3 yr Ϫ1 with ticks every 10 Ϫ3 yr Ϫ1 . All pdfs have been normalized to unity by dividing R by its zeroth moment, r 0 (⍀ i , ⍀ f ), which represents the mass fraction (here also the volume fraction) of the fluid elements that ventilate into the ocean through surface patch ⍀ i and out of the ocean through surface patch ⍀ f . In other words, what is shown are the ⍀ i → ⍀ f conditional residence-time pdfs. The numerical values of r 0 (⍀ i , ⍀ f ) are indicated in the figure in percent. (When r 0 is less than 0.1%, numerical values have been omitted for clarity.)

b. Model-derived residence-time distributions
To calculate R we first computed G and then diagnosed its flux through the surface. While in the continuum limit only the diffusive part of the surface flux, J · n, is nonzero by virtue of the boundary conditions on G, for the model we need to calculate the flux into the top 50-m-thick layer. Thus, numerically, J · n represents the vertical component of J through the 50-m-depth surface of both diffusive and advective fluxes. Physically, we can think of the fluxes through the 50-m surface as the fluxes through the base of the mixed layer. Figure 7 shows R if () ϵ R (, ⍀ i , ⍀ f ) for all combinations of entry and exit patches. The pdfs R if () can then be summed to obtain the pdfs of residence times from patch ⍀ i to anywhere on the entire surface ⍀ (rightmost column in Fig. 4), from anywhere on ⍀ to patch ⍀ f (bottom row in Fig. 4), or from anywhere on ⍀ to anywhere on ⍀ (bottom right-hand corner in Fig. 4; also Fig. 2). The general smoothness of the pdfs is a result of the spatial integration over our large surface patches and the time-averaged transport of the model.
If patches do not overlap or share a common edge, R if () is seen to go to zero, as → 0, as expected. Along the diagonal (identical in and out patches), and for the case where either the entry or exit patch consists of the entire surface, R if () approaches large values as → 0, with pairs (⍀ i , ⍀ f ) sharing an edge being intermediate in their short-behavior. The decomposition of R (, ⍀, ⍀) into the ⍀ i → ⍀ f contributions in Fig. 7 shows that the tendency of R (, ⍀, ⍀) to diverge for small comes from the short-circuit case of identical in and out patches (the diagonal in Fig. 7), consistent with the analysis in appendix A.
For several combinations of entrance and exit patches, R if () exhibits multimodality. Bimodality is particularly evident for fluid entering over the Southern Ocean (the ⍀ 5 row). We attribute the fast mode at a few decades to the relatively rapid transport of Antarctic Intermediate Water (AAIW). The second dominant mode at ϳ1000 yr and longer is consistent with the longer pathways from the abyss back to the surface, together with the fact that diapycnal diffusivities are weak. The model's ability to form distinct AAIW and Antarctic Bottom Water (AABW) was demonstrated by Primeau (2005). Also tabulated in Fig. 7 are the mass fractions of the ocean, r 0 (⍀ i , ⍀ f ), corresponding to the pdfs R ij (). Figure 7 displays pronounced asymmetry in that the ⍀ i → ⍀ f volume is very different from the ⍀ f → ⍀ i volume, which underscores the importance of advective transport. Purely diffusive transport would be self-adjoint with R if ϭ R fi .
As expected, the ocean's dominant water-mass formation regions are the high-latitude Atlantic (⍀ 1 ) and, especially for this model at least, the Southern Ocean (⍀ 5 ). More than 58% of the ocean's volume is ventilated into the Southern Ocean (⍀ 5 ) and approximately 40% of the ocean's volume is first exposed to the atmosphere in the Southern Ocean. A very large fraction of 28% [r 0 (⍀ 5 , ⍀ 5 )] is both ventilated in and out of ⍀ 5 , while only a few percent [r 0 (⍀ 1 , ⍀ 1 )] are ventilated both in and out of ⍀ 1 . Approximately 21% of the ocean's volume is formed in the high-latitude North Atlantic, and nearly 42% of this water (8.8% of the total volume) is first reexposed to the atmosphere in the Southern Ocean [r 0 (⍀ 1 , ⍀ 5 )]. The fraction doing the reverse, r 0 (⍀ 5 , ⍀ 1 ), is negligible, consistent with a global overturning circulation in which Atlantic surface waters flow northward. An important advantage of our diagnostics is that they account for both advective and diffusive transport. Conventional diagnostics, such as the meridional overturning streamfunction, ignore the important diffusive part of transport. Figure 8 shows the mean ⍀ i → ⍀ f residence time ϭ r 1 /r 0 in years. For most regions, the minimum value of lies along the diagonal (same in and out patches). A few off-diagonal patch pairs have shorter mean residence times but, with the exception of pathways involving the Southern Ocean (⍀ 5 ), they are insignificant in terms of the mass fraction of the ocean that they ventilate. Region ⍀ 5 is exceptional because of the pronounced bimodality of the ⍀ 5 → ⍀ f residence-time pdfs. This bimodality implies that the mean residencetime averages over two transport processes with widely separated modes (see Fig. 7). The relatively short residence times of the ⍀ 5 → (⍀ 4 , ⍀ 3 , ⍀ 2 ) pathways are manifestations of the strong first modes of the corresponding residence-time pdfs, which are dominated by the transport of AAIW to the Atlantic surface where it is first exposed to the atmosphere. The mean residence time of the ⍀ 5 → ⍀ 5 pathway is dominated by the much broader second mode of R 55 () located at տ 1000 yr. This second mode is the result of AABW filling up the abyssal ocean with slow transport back to the surface. In fact, the mean residence time of the ⍀ 5 → ⍀ 5 pathway is longer than for any other pathway with ⍀ 5 as the exit patch. One also notices that the ⍀ 5 → ⍀ 5 mean residence time of ϳ1300 yr is about an order of magnitude longer than any other mean residence times for identical entry and exit patches. This is consistent with the large 28% volume fraction that is ventilated into and out of ⍀ 5 , which implies that a significant number of fluid elements ventilated into the interior on ⍀ 5 journey far from ⍀ 5 before returning. Generally, the larger the distance between two patches, the larger , as expected.
The residence-time pdf would be a delta function with zero variance for a conveyor belt without any diffusion. The standard deviation, ϵ ͌ 2 Ϫ 2 , of the residence time is, therefore, a bulk measure of the amount of dispersion and mixing that occurs en route. Figure 9 shows a scatterplot of versus for all 121 combinations of in and out patches, which reveals a surprisingly compact relationship between and . The saturation of at ϳ1200 yr for the largest values may be attributed to the fact that the ocean has a finite volume and, hence, only a finite ability to disperse and mix fluid elements. Appendix B explains the saturation of at ϳ1200 yr in terms of the magnitude of the smallest eigenvalue 0 ϭ 1/797 yr Ϫ1 of the transport operator.

Discussion
To unambiguously characterize the rate of ventilation for any specified entry and exit patches, the full distribution as a function of residence time is needed. If one insists on summarizing the ventilation rate for overlapping entry and exit patches with a single number, some conditioning on the integration of must be imposed to eliminate the dominance of fluid with infinitesimal residence time. Possibilities include considering only fluid that will be transported at least a distance d 0 from the surface region of overlap or specifying a minimum residence time, 0 . For example, one might choose 0 ϭ 1 yr in an attempt to exclude that part of the ventilation flux due to waters that enter the seasonal thermocline but are reentrained into the mixed layer within the year before they reach the permanent thermocline.
However, merely cutting off the residence-time distribution below 0 does not cleanly separate seasonaland permanent-thermocline waters in the presence of diffusive mixing. Eddy-driven fluxes can transport water into and out of the permanent thermocline over a broad range of time scales, including times shorter than a year. In this study, we used an annually averaged transport operator that, as a result of this averaging, produces near-instantaneous vertical homogenization above the annual maximum mixed layer depth. Therefore, the ventilation rates computed with our steady ocean model should be interpreted as the rates with which waters enter or exit the permanent thermocline. The corresponding cumulative rate distributions for ventilation in or out of the permanent thermocline, ⌽(), show no plateaus or other natural breaks as a function of , implying sensitivity to 0 regardless of its value.
The residence-time pdf R and ventilation-rate distribution are flow diagnostics that can be applied to any type of model, from a plug-flow ventilated reservoir to a full OGCM and, once computed, do not need the model for their interpretation. This is in contrast to a summary of the domain-integrated vertical transport expressed as a single upwelling velocity and a single vertical diffusivity [e.g., the "recipes" of Munk (1966) and Munk and Wunsch (1998)], or as a single boxmodel-derived ventilation time scale or formation rate, all of which need the model to be interpreted. In this sense, R and , together with the boundary propagator on which they are based, are not only tracer-independent, but also model-independent transport diagnostics.
Notwithstanding the singularity at short residence times, there is another compelling reason why the ventilation rate is best described in terms of a distribution function. Transient tracers such as CFCs have been used for estimating residence-time-independent watermass formation rates Q from the formula where I is the oceanic inventory of a trasient tracer with surface concentration (t) (e.g., Broecker et al. 1998;Orsi et al. 1999). However, (18) only applies if the flow is purely advective so that ⌽ ↓ ϭ Q is independent of residence time . In general, eddy diffusive fluxes and mixing associated with entrainment give ⌽ ↓ a strong dependence on so that the tracer inventory must be obtained from the convolution Deconvolving (19) for ⌽ ↓ () by applying (18) is inappropriate for advective-diffusive flow because ⌽ ↓ is then not a constant that can simply be pulled out of the integral. If one does use (18) anyway, one does not Thus, Q is not a pure flow diagnostic but depends on the particular tracer employed. For example, Broecker et al. (1998Broecker et al. ( , 1999 used box models to obtain deepwater formation rates of 5 Sv from CFC-11 inventories and of 15 Sv from 14 C distributions. Based on this result Broecker et al. (1999) went on to speculate that there has been a major slowdown of deep-water formation during the twentieth century. In a subsequent paper, Orsi et al. (2001) showed that the disparate estimates of 5 and 15 Sv can be reconciled without invoking a major slowdown in deep-water formation rates if the CFC and 14 C data are interpreted using a two-box model instead of the one-box model used by Broecker et al. (1998).
Our work suggests that the agreement or disagreement between various estimates of deep-water formation may well be an artifact of inadequate representation of eddy diffusive fluxes. Quasi-random mesoscale eddies are believed to be very important in the watermass formation process for both open ocean and coastal settings (Marshall 1997;Marshall and Schott 1999;Pringle 2001). Such processes are reasonably modeled as continuum diffusion: their outright neglect cannot be justified.
Although our discussion has focused primarily on using transient tracers to estimate ventilation rates, it is important to point out the connection of our work to the problem of estimating rates of water-mass formation and subduction using the so-called thermodynamic and kinematic methods (e.g., Walin 1982;Marshall et al. 1993;Speer and Tziperman 1992;Marshall et al. 1999). In these studies, the circulation is assumed to be steady on time scales longer than a year, and the buoyancy flux through the mouth of out-cropping isopycnals is estimated either (i) by using surface fluxes of heat and moisture to estimate a rate of water-mass transformation whose convergence drives subduction (thermodynamic method) or (ii) by directly combining windderived Ekman pumping rates with the geostrophic flow across the sloping base of the mixed layer to obtain a subduction rate [kinematic method; Marshall et al. (1993)]. Ventilation rates based on the kinematic method neglect the important contribution of eddydriven diffusive fluxes in ventilating the ocean interior (e.g., Marshall et al. 1993) and, to the extent that the other methods include diffusive fluxes, they focus on the net flux of buoyancy and not on determining a ventilation rate that can be used to estimate the invasion of an arbitrary tracer into the ocean.

Conclusions
The ocean is a vast reservoir of atmospheric constituents that enter the ocean with a characteristic global pattern of fluxes, remain dissolved in the oceans for times ranging from seconds to millennia, and are then reventilated to the atmosphere with another characteristic pattern of fluxes. This paper presents a new approach for thinking about oceanic ventilation rates. We have developed ventilation-rate diagnostics derived from the boundary propagator G and its surface fluxes. Because G is a type of Green function that is an integral representation of the fluid's transport operator, our diagnostics are tracer independent and quantify the surface-to-surface transport of oceanic flow.
We have introduced ⌽ ↓ (⍀ i , t i ; ⍀ f , ), the flux at time t i through a surface patch ⍀ i of newly ventilated water that is destined to be reexposed to the atmosphere at ⍀ f after having resided in the interior for a time or longer. Symmetrically, ⌽ ↓ (⍀ f , t f ; ⍀ i , ) is the outward flux through ⍀ f at time t f of water that was formed at ⍀ i and resided in the interior for a time or longer. In terms of , ⌽ ↓ and ⌽ ↑ are cumulative ventilation-rate distributions. The corresponding ventilation-rate density distributions, ↓ ϭ Ϫd⌽ ↓ /d and ↑ ϭ Ϫd⌽ ↑ /d, give the ventilation rate per unit residence time. The distributions and R are intimately related through a spectral generalization of the familiar result that the residence time in a reservoir is the ratio of the mass in the reservoir to the flux through the reservoir.
Eddy diffusive fluxes impart a short-singularity to R () and ⌽() if the surface region where fluid enters the ocean interior overlaps with the region where fluid is reexposed to the atmosphere. For overlapping entry and exit regions, R and ⌽ diverge like Ϫ1/2 as → 0. While the pdf R is always integrable and properly normalized, the singularity of ⌽ renders the ventilation-rate density distribution ϭ Ϫd⌽/d unintegrable. This means that, for overlapping entry and exit regions the ventilation rate regardless of the ventilated fluid's residence time in the interior, ͐ ϱ 0 ()d ϭ ⌽(0), is infinite and, therefore, cannot meaningfully be estimated.
We have used an off-line OGCM to illustrate the rich detail that our diagnostics provide on the way that the ocean communicates with the surface. When ⍀ i and ⍀ f are the entire ocean surface, ⌽() represents the global inventory of water binned according to age. This global ⌽() and the corresponding residence-time pdf, R (), have a remarkably simple structure with near-powerlaw scaling from ϳ10 to ϳ600 yr. The complex nature of the model's surface-to-surface transport is revealed only when ⌽ or R are partitioned spatially according to where fluid elements ventilate into and out of the ocean interior.
We have quantified which geographical locations play an important role in ventilating the ocean interior. The temporal norm of R gives the fraction of the ocean volume involved in connecting any pair of specified entry and exit patches. Maps of the flux of newly ventilated water into the interior identify the subtropical gyres, the high-latitude North Atlantic, and the Southern Ocean as the major water-mass formation regions. Maps of the exposure flux of old water identify equatorial bands and the Southern Ocean as the major upwelling regions. The Southern Ocean stands out, at least in our model, as playing a dominant role both in the formation of newly ventilated water and in the reexposure of older waters. The residence-time distributions for waters formed in the Southern Ocean, and first reexposed to the atmosphere in the Atlantic, are bimodal. A mode at a few decades is associated with relatively fast transport of Antarctic Intermediate Water, and a second broader mode at ϳ1000 yr is associated with longer transport pathways to the abyss and back to the surface.
Comparing our model results with published ventilation rates obtained from observational data is difficult because the latter are not usually tracer independent, and their implied conditioning of the ocean's true ventilation-rate distribution is not explicitly controlled. It will be imperative to explore the constraints imposed by oceanic tracer data on the ventilation-rate distributions and to reinterpret published ventilation rates in terms of these distribution.
While we have not examined the interior pathways of surface-to-surface transport, this is an interesting area of future research given the current interest in understanding the relative roles of surface forcing and crossisopycnal eddy mixing in driving the meridional overturning circulation, and the relative dearth of obser-vations on how deep water returns to the surface (Toggweiler and Samuels 1998;Wunsch and Ferrari 2004). Future work will focus in detail on the three-dimensional surface-to-surface paths in terms of the extended pdf of residence time, ventilation locations, and interior transit location. Investigating the sensitivity of our model-derived ventilation rates to changes in surface momentum forcing, buoyancy forcing, and subgrid-scale mixing was beyond the scope of the present study, but is an important avenue of future research. Because our diagnostics provide an integrated measure of transport to and from the sea surface, they are also well suited to assessing the impact of changes in different physical processes on modeled ocean biogeochemical cycles. Furthermore, we plan to exploit our diagnostics to quantify how water-mass formation rates change with the seasons (e.g., Williams et al. 1995) and from year to year as the result of climate change.
as → 0, where A i is the area of ⍀ i . While the prefactors of (A4) are rough estimates, it is clear that for ⍀ i ϭ ⍀ f , the distribution ↓ must diverge like Ϫ3/2 as → 0.
If ⍀ f and ⍀ i are not identical, but merely overlap, the flux out of ⍀ f is not approximately the flux out of the entire ocean, but reduced from that by a factor of A overlap /A i . Otherwise, the argument above for a Ϫ3/2 divergence follows through, but with a prefactor reduced from that in (A4).
If ⍀ f and ⍀ i merely share an edge, then the flux through ⍀ f is effectively reduced from the flux through the entire ocean by A diff /A i , where A diff is a strip of ⍀ f along its edge with ⍀ i that has a width on the order of ͌ and a length on the order of ͌A i . (For small most of the flux is through ⍀ i except for some flux through a diffusively accessible ribbon along the edge with ⍀ f .) Thus, to a first approximation, we can multiply the time derivative of the volume-integrated G in (A3) by (͌A i ͌)/A i to capture the effective fraction of the global flux. As → 0, the diverging flux is therefore compensated by a shrinking area of ⍀ f that is diffusively accessible so that, in this case, ↓ → Ϫ1 .

Asymptotic Standard Deviation of Residence Time
One of the notable features in Fig. 9 is the saturation of the standard deviation at approximately 1200 yr for mean residence times Ͼ 1500 yr. This appendix provides an explanation for this behavior.
Suppose one were to randomly select a group of N paths from ⍀ i to ⍀ f with residence times ( 1 , 2 , . . . , N ) and compute the sample mean, , and the sample standard deviation, . We will show that, as N → ϱ, the expected value of conditioned on the sample mean taking on the particular value, ϭ , satisfies where is the smallest eigvenvalue of the transport operator.
The joint pdf describing the sample residence times is given by f͑ 1 , 2 , · · ·, N ͒ ϭ R ͑ 1 ͒R ͑ 2 ͒ · · · R ͑ N ͒, ͑B2͒ where R () is shorthand notation for R (, ⍀ f , ⍀ i ). The conditional pdf describing the sample given that ϭ is given by f evaluated along the N Ϫ 1 dimensional plane ( 1 ϩ 2 ϩ · · · ϩ N ) ϭ N. In the limit that → ϱ, we need only consider the tail of R when evaluating f because the fraction of the total "area" on the hyperplane that is in the strips 0 Յ i Ͻ 0 for i ϭ 1, 2, . . . , N approaches zero. In other words, for large it is increasingly unlikely that any of the i in the sample will be less than 0 .
The tail of (6) and (7)] is determined by the tail of ⌽(), and we showed ⌽() to be the mass-weighted domain integral of G. Because the ocean has finite volume, the transport operator has a discrete eigenmode expansion, guaranteeing that G, and hence ⌽, decay asymptotically like exp(Ϫ 0 ), where 0 is the smallest eigenvalue of the transport operator. The tail of R () is thus given by R ͑͒ ϳ 0 2 exp͑Ϫ 0 ͒ as → ϱ, ͑B3͒ and the conditional joint pdf describing our sample satisfies lim →ϱ f͑ 1 , 2 , . . . , N | ϭ ͒ ϭ T͑ 1 ͒T͑ 2 ͒ · · · T͑ N ͒,

͑B4͒
where T͑͒ ϭ 0 2 exp͑Ϫ 0 ͒ ͑B5͒ is the normalized tail distribution [expression (B3) happens to be already normalized]. Heuristically, for particles that reside in the ocean for a very long time, the ocean begins to look like a well-mixed box. As N → ϱ, will approach the standard deviation of T (). A straightforward calculation shows that the standard deviation of T () is ͌2/ 0 . For our steady-state transport model, 1/ 0 ϭ 797 yr (obtained by applying the inverse power method to the transport operator.) Thus, the expected asymptote for the standard deviation is 1127 yr, consistent with the large-points in Fig. 9.