Characterizing Transport between the Surface Mixed Layer and the Ocean Interior with a Forward and Adjoint Global Ocean Transport Model

The theory of first-passage time distribution functions and its extension to last-passage time distribution functions are applied to the problem of tracking the movement of water masses to and from the surface mixed layer in a global ocean general circulation model. The first-passage time distribution function is used to determine in a probabilistic sense when and where a fluid element will make its first contact with the surface as a function of its position in the ocean interior. The last-passage time distribution is used to determine when and where a fluid element made its last contact with the surface. A computationally efficient method is presented for recursively computing the first few moments of the first- and last-passage time distributions by directly inverting the forward and adjoint transport operator. This approach allows integrated transport information to be obtained directly from the differential form of the transport operator without the need to perform lengthy multitracer time integration of the transport equations. The method, which relies on the stationarity of the transport operator, is applied to the time-averaged transport operator obtained from a three-dimensional global ocean simulation performed with an OGCM. With this approach, the author (i) computes surface maps showing the fraction of the total ocean volume per unit area that ventilates at each point on the surface of the ocean, (ii) partitions interior water masses based on their formation region at the surface, and (iii) computes the three-dimensional spatial distribution of the mean and standard deviation of the age distribution of water.


Introduction
The interaction of the ocean with the rest of the climate system happens primarily at the sea surface through air-sea fluxes. These air-sea fluxes imprint the current physical and chemical state of the atmosphere on water parcels. When water parcels are transported below the surface, they become shielded from the atmosphere until they resurface and communicate past physical and chemical climate conditions to the atmosphere. The delay between successive visits to the mixed layer imparts to the earth system an important long-term memory. Understanding the duration of this memory is more complicated than the simple "conveyor belt" picture of the global ocean circulation, which suggests that fluid parcels are transported below the surface at one end of the conveyor belt, carried steadily along, and fatefully popped up at the other end. In reality fluid elements are transported to and from the mixed layer by advective, mixing, and diffusive processes that lead to a continuous distribution of times between successive visits to the surface mixed layer. Since fluid element properties are reset by air-sea interactions during each visit to the surface mixed layer, the "memory" of a fluid element does not precede its last contact with the surface layer, nor does it persist beyond its next contact with the surface layer. Ocean transport is therefore described by the distribution of times for fluid elements to have had their last contact with the atmosphere as well as the distribution of times when fluid elements will have their first (or next) contact with the atmosphere.
The analysis of geophysical flows in terms of distribution of times for fluid elements to be transported from a surface of interest to the interior of the fluid body is a growing subject. A useful review in the context of stratospheric air can be found in Waugh and Hall (2002). In an oceanographic context, Haine and Hall (2002) combined the familiar concepts of water mass and age into a generalized theory formulated in terms of the surface boundary condition propagator (boundary Green function) of the ocean's advectiondiffusion transport equation. The connection between the forward and adjoint transport equations and the distribution of times for fluid elements to be transported to and from a surface is discussed in Holzer and Hall (2000). Delhez and Deleersnijder (2002), formu-late a more general theory that allows the age of a tracer constituent to be tracked. The connection between the surface-to-interior transit time distribution and the different concepts of tracer age used by oceanographers is discussed in Hall and Haine (2002) and in Waugh et al. (2003).
In this paper we explore when, where, and how much fluid is exchanged between the surface mixed layer and the interior ocean in a global ocean general circulation model (OGCM) by computing first-and last-passage time distributions. 1 The quantities we compute are not only useful diagnostics for characterizing the ventilation of an OGCM but are properties of the ocean that are not directly observable but that are nevertheless of direct interest to oceanographers studying the exchange of CO 2 between the ocean and the atmosphere. The focus of this paper is not on the transport of any particular tracer, but rather on the physical transport of fluid elements caused by currents and mixing processes.
The plan of the paper is as follows. In section 2 we introduce the first-and last-passage time distributions and their connection to the residence time of a fluid parcel in the ocean interior. The transport problems that must be solved to compute the distributions and their moments are also given in this section. In section 3 we describe the offline ocean transport model and the dynamical OGCM. Section 5 presents surface maps that show the volume fraction of the total ocean volume that ventilates at each surface point of the ocean. Section 6 presents a watermass analysis for the model that was computed using the last-passage time distribution. Section 7 presents an age analysis as well as an analysis of the time for water in the interior to return to the ocean surface. Section 7a discusses the shape of typical last-passage time distributions. Section 8 discusses the results and concludes.

First-and last-passage time distributions
We begin by considering a fluid element located below the sea surface in the ocean interior, and denote by r , the total amount of time the fluid element will spend in the ocean interior-that is, the amount of time between its last contact with the sea surface and its next contact with the sea surface. If we consider a fluid el-ement in the ocean interior at time t, its residence time can be decomposed into two parts, where the first part, lp , is the time that has elapsed since the fluid element made its last contact with the sea surface and the second part, fp , is the time that will pass before the fluid element makes its next contact with the surface. We refer to fp as the first-passage time in accordance with the terminology used in the literature on stochastic processes (e.g., Stratonovich 1963). By analogy, we will refer to lp as the last-passage time. The last-passage time is also refered to as the transit-time or age in the oceanographic literature (e.g., Haine and Hall 2002;Wunsch 2002). Because of diffusive and mixing processes, the fluid elements in a parcel are transported to and from the surface via many pathways, and a distribution is needed to describe when and where the fluid elements in a water parcel make their last and first contact with the sea surface.
To make these ideas precise, we will denote the ocean surface by ⍀ s , a point on the surface by x s , and a point in the ocean interior by x (Fig. 1). The joint probability density that a fluid element was at x s on ⍀ s at time t Ϫ lp , given it is currently at interior position x at time t, is denoted by L(x s , lp | x, t). The joint pdf that a fluid element at position x at time t will make its first contact with the surface at time t ϩ fp and and position x s , will be denoted by F(x s , fp | x, t).
In the following, we will focus on the situation where the circulation is assumed to be steady. In other words, we will assume that the transport by time-dependent eddies can be represented by a stationary eddy diffusion tensor. For the stationary case, the joint distribution for (x s , lp ) and (x s , fp ) are independent of t and are thus conditioned only on the interior position x.

a. Last-passage time problem formulation for a stationary flow
As shown in Holzer and Hall (2000) and in Haine and Hall (2002), the joint distribution, L(x s , lp | x), for the last-passage time, lp , and last surface-contact position, x s , can be expressed in terms of the propagator of surface boundary conditions L( subject to the surface boundary condition, In (2) T is the steady advection-diffusion transport operator where U is the fluid velocity and K is the diffusion tensor, and in (3) Physically, Eq.
(2), together with the surface boundary condition (3), describes the time evolution of tracer-tagged fluid elements as they are transported into the ocean interior. The fluid elements are tagged with tracer when they make contact with the surface at position x s and time t s and subsequently untagged when they make their next contact with the surface.
b. First-passage time problem formulation for a stationary flow Holzer and Hall (2000) also show that the joint probability density, F(x s , fp | x), for the first-contact position x s and first-passage time fp to ⍀ can be expressed in terms of the propagator for surface boundary con- . Note that Eq. (5) is to be integrated with respect to t backwards in time starting from t ϭ t s .

c. Moments
The first-and last-passage time distribution functions F and L contain all information about transport to and from the sea surface in time-integrated form. Even for the case of stationary flow, the joint distributions, L(x s , lp | x) and F(x s , fp | x), can contain an overwhelming amount of information. For example, in a discrete ocean model with n x ϫ n y surface points and n z vertical levels, the L and F each contain n 2 x ϫ n 2 y ϫ n z continu-ous functions of time. To digest this vast amount of information, it is useful to summarize the joint distributions in terms of various moments and integrated quantities that can be interpreted physically. Most of these quantities will be introduced and interpreted in the following section. Here we define some of the more fundamental quantities that are used throughout the rest of the paper. The marginal distributions obtained by integrating out x s , yield the probability density for the first-and lastpassage time to the surface ocean ⍀ s taken as a whole. Useful summary statistics for L and F are the mean and standard deviation. For the last-passage time distribution, the mean (⌫) and standard deviation (⌬) are given by Similarly, the mean and standard deviation of the firstpassage time distribution The mean last-passge time, ⌫(x), is often refered to as the mean age or ideal age (Holzer and Hall 2000;Haine and Hall 2002;Khatiwala et al. 2001). The age terminology comes from thinking of a fluid parcel as having been born when it made its last contact with the surface. By analogy, if we think of the first-contact with the surface as the death of the fluid parcel, then ⌫ (x) is the remaining life expectancy of the fluid parcel at x. The age gives a characteristic time scale for the transport from the sea surface to the interior positon x. Because first-and last-passage time distributions have finite width, a single time scale is inadequate for describing the transit time to or from the sea surface (Wunsch 2002;Waugh et al. 2003). Because of this, the width of the distribution is a useful complement to the mean age. The standard deviations ⌬ and (⌬ ) measure the width of the last-and first-passage time distributions. They can be interpreted as a measure of the amount of mixing a fluid parcel experiences en route from and to the sea surface (Delhez and Deleersnijder 2002;Waugh et al. 2003). It is easy to show (see, e.g., Stratonovich 1963) that by multiplying Eqs. (2) and (5) by n lp and n fp , respectively, and integrating by parts with respect to lp and fp the moments of the distributions can be obtained recursively by solving time-independent boundary value problems. For example, ⌫ is given by the solution to subject to ⌫ ϭ 0 at the surface, and ⌬ can be obtained by solving for the second moment of the distribution (m 2 ) using the following equation: subject to m 2 ϭ 0 at the surface and setting ⌬ ϭ ͌m 2 Ϫ ⌫ 2 . To compute ⌫ and ⌬ , the adjoint transport operator T † is used instead of T.

Ocean transport model
In this section we apply our theory to the timeaveraged transport operator from a global ocean general circulation model. We use a dynamical model described below to produce the three-dimensional flow and mixing fields. These fields are then used as input to an offline transport model. The transport model is coded in Matlab using a sparse matrix formulation that allows us to easily obtain the adjoint model and to efficiently invert the forward and adjoint transport operators to compute the moments of the distribution. To solve the stationary transport problems, we used MUMPS (Amestoy et al. 2001) to factor the transport operators. The full distributions presented in section 7a were computed by writing the discretized solution of Eqs.
(2)-(3) in the form of a matrix exponential and using the Expokit software package (Sidje 1998) to evaluate the solution at different times. The advantage of the matrix exponential approach is that it is well suited for stiff problems where large diffusivities enter the transport operator in regions of convective overturning. The method also handles the convergence of the mesh near the poles without requiring filtering.

Ocean model configuration
We use a version of the ocean component of the Canadian Centre for Climate Modeling and Analysis climate model (NCOM). The model has 29 levels in the vertical ranging in thickness from 50 m near the surface down to 300 m near the bottom. The horizontal resolution has 48 grid points meridionally and 96 grid points zonally for an approximate resolution of 3.75°ϫ 3.75°. The model uses the KPP (Large et al. 1994) vertical mixing scheme and the GM (Gent and McWilliams 1990) isopycnal mixing scheme. Tracers are advected using a second-order centered difference scheme. The dynamical model is forced by prescribed monthly freshwater and heat fluxes (obtained from the atmospheric component of the climate model) together with a restoring to observed surface temperature with a 30-day time scale and to observed surface salinity with a 180day time scale.
The model was spun up for over 8000 years. At the end of the spinup period, the flow and diffusion tensor fields were averaged over a period of 5 yr. The timeaveraged fields were then used to construct the timeaveraged advection-diffusion transport operator, T. The time-averaged operator was used for all computations presented in this report.

Partition of surface ocean
We will, on several occasions, present diagnostics that summarize the transport to and from a set of discrete nonoverlapping surface regions. It is convenient to introduce these regions now. Figure 2 shows how the surface of the ocean is partitioned into 11 discrete regions that together cover the whole surface of the ocean: This partition is similar to the one used to define regions of ocean-atmosphere carbon dioxide exchange in the TransCom project (Gurney et al. 2000).

Characterizing surface ocean points based on their importance for ventilating the interior of the ocean
The sites at the sea surface where "new" water is formed and where "old" water is exposed to the atmosphere are important characteristics of the ocean ventilation process. Knowing how much water is ventilated in any surface region of the ocean can contribute to our understanding of the cycling of chemical constituents between the ocean and atmosphere. For example, the efficiency of the "biological pump," the sequestration of atmospheric CO 2 into the ocean interior by sinking organic carbon particles, depends on how much of the upwelled (preformed) nutrients are utilized by organisms before being returned to the ocean interior through the formation of new water (Sigman and Boyle 2000). To evaluate this efficiency, one must compute the globally averaged surface concentration of preformed nutrients weighted by the volume of the ocean that is ventilated through each surface patch. As far as we know, evaluating the efficiency of the biological pump in this way has never been done for a global OGCM. The likely reason is that it is difficult to determine what volume of the ocean is ventilated at each point on the ocean surface. In this section we present an efficient method for determining how much water is ventilated at each point on the surface of the ocean. We will also apply the method to the global OGCM described in section 3a and present maps that show where new water is formed and where old water is first exposed to the atmosphere.
It is useful to introduce the thickness H new (x s ), defined such that the volume of the ocean that gets formed through an infinitesimal surface element dx s centered at x s , can be expressed as , is a volume density (volume per unit area). It can be obtained by integrating out lp and x in the last-passage time distribution L(x s , lp |x): In the above, the function H new (x s ) is an inventory of the entire ocean volume based on where each fluid element made its last contact with the surface. We can also form an inventory of the ocean based on where water will be first exposed to the atmosphere. This inventory can be expressed in terms of the volume density H old (x s ) by integrating out fp and x in the firstpassage time distribution: It turns out that both H old and H new also have interpretations in terms of fluxes: H old is also the outward flux of an ideal age tracer ⌫(x) (mean last-passage time) through the surface at x s , Analogously, H new is the outward flux through x s in the time-reversed adjoint flow of the mean first-passage time ⌫ (x), A formal derivation of these results follows readily from the reciprocity relationships for the forward and adjoint transport equations described in Holzer and Hall (2000). We hope to develop these ideas further in a follow-up paper, but an intuitive understanding of the flux relationships can be obtained from a demographic analogy in which we think of a fluid particle's last-contact with the sea surface as its birth and of its firstcontact with the sea surface as its death. In this analogy, lp is the age of the particle and fp is its remaining lifetime. The function H new (x s )dx s is the size of the population of particles in the ocean interior that were born in the surface patch dx s centered on x s , and H old (x s )dx s is the size of the population of particles in the ocean interior that will die in the patch centered on x s . With the demographic analogy in mind, (16) determines the population size H new (x s )dx s from a direct census over the whole ocean volume, while the flux relationship (18) takes advantage of the stationarity of the age distribution to determine the same population size from the birthrate at x s weighted by the life expectancy of the particles born there. Similarly, the flux relationship (17) estimates the population size, H old (x s )dx s , from the deathrate at x s weighted by the age of the particles that die at x s .
Apart from having an intuitive interpretation in terms of the ventilation of new and old water, the flux relationships (18) and (17) have an important computational advantage over the direct census approach of Eqs. (15) and (16); they allow H new (x s ) to be computed using a single "ideal age" tracer (⌫) and H old (x s ) to be computed using a single "remaining life expectancy" tracer (⌫ ). Obtaining the same results by first solving Eq. (2) subject to the boundary condition (3) and then using (15) would require a separate time-dependent tracer simulation for each surface grid point-a prohibitively large computation.
We applied the above method to the transport operator of the global OGCM. The steps used to obtain H old were to first compute the mean of the last-passage time distribution, ⌫, by solving T ⌫ ϭ 1 subject to ⌫ ϭ 0 at the surface, and then computed the flux of ⌫ through the surface to obtain H old according to Eq. (17). Similarly, to obtain H new we first compute the mean of the first-passage time distribution, ⌫ , by solving T † ⌫ ϭ 1 subject to ⌫ ϭ 0 at the surface, and then compute the flux of ⌫ through the surface using the adjoint flow field, according to Eq. (18). See section 7 for a description of the ⌫(x) and ⌬(x) fields.
The new and old water thickness results are shown in Fig. 3. The inventory of water in the interior ocean binned according to where it made its last contact with the surface is shown in the top panel and the inventory of the ocean volume based on where the water will make its first contact with the surface is shown in the bottom panel.
The formation of new water is localized to a few small regions in high latitudes with the biggest contributions coming from locations in the North Atlantic Ocean and the Weddell Sea. The contributions from the downward Ekman pumping regions in the subtropical gyres are relatively small. By comparing the upper and lower panels of Fig. 3, we see that the regions where older water is first exposed to the atmosphere are broader than the formation regions but are never- theless nonuniform as is assumed in the Stommel-Arons theory of the abyssal circulation (Stommel and Arons 1960). Upwelling regions in the Southern Ocean, the eastern equatorial Pacific Ocean and at the midlatitude western boundaries are the dominant regions for older water to make its first contact with the surface. It is important to point out that our calculation includes both the effects of vertical advection and eddy diffusion in the ventilation process. This is in contrast to the work of Qiu and Huang (1995) that studied only the advective component of the ventilation process.
In Table 1 we present a summary of the fraction of the total ocean volume that ventilates into and out of the ocean interior in each surface region shown in Fig.  2. The dominant water formation regions are the North Atlantic (⍀ 1 ), responsible for 21% of the total ocean volume, and the Southern Ocean (⍀ 5 ), responsible for 58% of the total volume. The dominant region where older interior water makes contact with the surface for the first time is the Southern Ocean (⍀ 5 ), with 40% of the total volume. The North Pacific, region ⍀ 11 , is also important because of its large surface area; 15% of the total ocean volume makes its first contact with the surface in region ⍀ 11 . The upwelling region of the western equatorial Pacific (⍀ 10 ) is also an important region for older water to become exposed to the atmosphere. Approximately 11% of the total ocean volume is ventilated in a region that accounts for only 6.4% of the total ocean surface area.

Characterizing interior water parcels based on
where they were formed at the surface (watermass analysis) Determining the surface source regions for water in the interior of the ocean is related to the problem of watermass analysis. Haine and Hall (2002) make a precise connection between the concept of water mass and the Green function for surface boundary condition by interpreting the boundary Green function as a joint distribution of composition and time of last contact with the surface. In this formalism, the traditional concept of water mass based only on composition is recovered by integrating out lp in L(x s , lp |x) and further integrating the resulting marginal distribution over a surface patch, ͑19͒ The quantity f ⍀ i (x) is the volume mixing ratio at x of water that was ventilated into the ocean interior from the surface region ⍀ i . If the surface of the ocean is partitioned into a finite set of nonoverlapping surface regions,   Fig. 2. The first column shows the fraction of the total surface area covered by each region. The second column shows how much new water (f ⍀ i ) is formed in each region (last contact). The third column shows how much old water (f ⍀ i ) is exposed to the atmosphere in each region (first contact).

Region
Surface area fraction (%) water parcels can be decomposed into a finite number of watermass types such that The n volume mixing ratios, f ⍀ i (x), account for all of the water in a parcel. Maps and sections of f ⍀ i (x) are useful for analyzing the interior distribution of water that was ventilated into (formed in) region ⍀ i . A related set of quantities, for i ϭ 1, 2, · · · , n, partition a water parcel at x based on where the fluid elements contained in it will make their first contact with the surface. The quantities, f ⍀ i (x) for i ϭ 1, 2, · · · , n, also sum to unity and account for all of the water in the parcel. Maps and sections of f ⍀ i (x) tell us what fraction of the water will ultimately become exposed to the atmosphere in region ⍀ i . The volume mixing ratio f ⍀ i can be obtained directly by solving the forward transport problem The adjoint volume mixing ratio f ⍀ i is obtained from the analogous problem in which the adjoint transport operator, T † , is used instead of T.
We have applied the above diagnostic to the OGCM with a partition of the surface into the n ϭ 11 regions shown in Fig. 2. We computed f ⍀ i (x) by solving a discretized version Eq. (23) for each ⍀ i . In the discretized problem, the surface boundary condition was applied to the uppermost model level so that the effective ⍀ s surface is at the bottom of the uppermost layer of the model, which corresponds to a depth level of 50 m.
In the following sections we describe meridional cross sections of f ⍀ i (x) taken from each ocean basin.
a. Atlantic Figure 4 shows depth-longitude cross sections of f ⍀ i (x) along 30°W in the Atlantic Ocean. The contours indicate the volume fraction of water that made its last contact in each surface region. We only show those surface regions that have a significant contribution to the volume mixing ratio somewhere along the section.
A broad tongue of water that made its last surface contact in region ⍀ 1 in the high-latitude North Atlantic extends southward as can be seen with the f ⍀ 1 ϭ 0.5 contour that extends up to 40°S at a depth of 1872.5 m (top-left panel). This water corresponds to North Atlantic Deep Water (NADW). The water above and below the NADW tongue originates predominantly in the surface of the Sourthern Ocean (region ⍀ 5 ). Source region ⍀ 5 (bottom-left panel) contributes a core of water that spreads northward in an intermediate depth range centered at a depth of approximately 550 m. This core of water corresponds to Antarctic Intermediate Water. Source region ⍀ 5 is also the dominant contributor to the deep waters of the Atlantic. More than 80% of the bottom water below 3500 m, corresponding to Antarctic Bottom Water, originates from region ⍀ 5 . (See also the discussion relating to 14 C/C observations in section 7.) The surface regions of the subtropical gyres in the Northern Hemisphere (region ⍀ 2 : Fig. 4, top-right panel) and Southern Hemisphere (region ⍀ 4 : middleright panel) are significant contributors to water in the top 1100 m. Other surface regions from the Pacific and Indian Oceans do not contribute significant amounts of water to the Atlantic Ocean.

b. Indian Ocean
The panels of Fig. 5 show meridional cross sections of f ⍀ i (x) at 90°E in the Indian Ocean. The contours indicate the volume fraction of water that made its last contact in each of the different surface regions. At this longitude, water that made its last surface contact in region ⍀ 1 of the North Atlantic (upper-left panel) has a core maximum of more than 25% near 30°S at a depth level 2350 m and spreads horizontally as indicated by the slab of water in the depth range between 1300 and 2800 m confined by the f ⍀ 1 ϭ 0.2 contour.
Water that made its last surface contact in the Southern Ocean (region ⍀ 5 ) contributes more than 80% by volume of the water in the bottom of the Indian Ocean (Fig. 5, upper-right panel). This water mixes upward, apparently filling the basin from below.
In the latitude band directly below region ⍀ 6 , more than 90% by volume of the water in the thermocline above 700 m made its last contact with the surface in region ⍀ 6 (middle-left panel). Contributions from region ⍀ 6 also extend northward, with more than 10% of the water above 1400 m and below 100 m having had its last surface contact in region ⍀ 6 . Contributions from region ⍀ 8 (bottom-left panel) is confined in the upper thermocline with the f ⍀ 8 ϭ 0.1 contour extending to a maximum depth of approximately 650 m in the northernmost part of the section.
Water from the surface of subtropical south pacific (region ⍀ 7 ) entering the Indian Ocean through the Indonesian Throughflow is also a significant contributor at 90°E. Source water from region ⍀ 7 has a core maximum of more than 15% by volume at a depth of 250 m near 10°S.

c. Pacific Ocean
The panels of Fig. 6 show meridional cross sections of f ⍀ i (x) at 160°W in the Pacific Ocean. The contours indicate the volume fraction of water that made its last contact in each of the different surface regions.
There is a core of water originating from region ⍀ 1 that enters the Pacific via the Southern Ocean (upperleft panel) as indicated by the f ⍀ 1 ϭ 0.2 contour in the southern part of the section. This water also makes a relatively uniform contribution of 15%-20% by volume of the water below the thermocline. In the thermocline, the contribution from region ⍀ 1 gradually decreases as one moves upward toward the surface.
Water that made its last surface contact in the Southern Ocean (region ⍀ 5 ) is the dominant contributor in the deep Pacific (upper-right panel). More than 70% by volume of the water below approximately 2000 m was last at the surface in region ⍀ 5 . Contributions from ⍀ 5 decrease rapidly as one moves upward into the thermocline. Water that made its last surface contact in regions ⍀ 7 and ⍀ 11 is confined to the thermocline (middle-left panel and bottom-right panel). There is also a significant cross-equatorial transport for water that originated in regions ⍀ 7 and ⍀ 11 . The f ⍀ 11 ϭ 0.05 contour for water that made its last surface contact in the subtropical South Pacific (region ⍀ 7 ) extends north- FIG. 4. Meridional cross section in the Atlantic at 30°W of f ⍀i , showing the probability that the fluid along the section had its last surface contact in regions ⍀ 1 , ⍀ 2 , ⍀ 3 , ⍀ 4 , or ⍀ 5 . Only the regions that contribute at least 10% of the volume fraction at some point in the section are plotted. The solid contour interval is 0.1, and the dotted contour interval is 0.05. ward of 25°N. Similarly, the 5% contour for water that made its last surface contact in the subtropical and subpolar North Pacific (region ⍀ 11 ) extends southward to approximately 25°S. Water that made its last surface contact in the equatorial Pacific (region ⍀ 9 : middleright panel, and region ⍀ 10 , not shown) is confined to a thin lens close to the surface with the f ⍀ 10 ϭ 0.1 contour reaching no deeper than 250 m.
There is also a core of water that made its last surface contact in the south Indian Ocean (region ⍀ 7 : Fig. 6, bottom-left panel). This core lies above the core of water originating from region ⍀ 1 in the North Atlantic and has a maximum of more than 15% by volume near 45°S at a depth of approximately 950 m. This water also spreads northward along isopycnal surfaces as indicated, for example, by the 5% contour that extends northward to approximately 10°N.

Characterizing the time scales for the transport to and from the surface (age analysis)
The concept of "age," the time elapsed since water had last contact with the atmosphere, is invaluable for FIG. 5. As in Fig. 4, but for a meridional cross section in the Indian Ocean at 90°E.

APRIL 2005 P R I M E A U
interpreting ventilation rates estimated from transient and radio tracers. Many studies have used oceanic tracer data to estimate the "age" of water parcels in the ocean interior (e.g., Jenkins 1998). Nevertheless, because of mixing and diffusion, there is no single time scale that best characterizes when a fluid parcel was last exposed to the atmosphere; instead there is a distribution of ages. The width, ⌬ of the distribution is, in fact, of crucial importance for analyzing age fields and transient tracer data in fluid systems with mixing (Delhez and Deleersnijder 2002;Waugh et al. 2003; see also Hall et al. 2004 for an important application to the problem of estimating the oceanic uptake of anthropogenic CO 2 ). While there have been several modeling studies that have studied the mean (⌫) of the age distribution as a way of characterizing the ventilation in GCMs (England 1995;Goodman 1998;England and Holloway 1998;Khatiwala et al. 2001), few studies, if any, have analyzed the width of the age distribution simulated in a global OGCM. We therefore present FIG. 6. As in Fig. 4, but for a meridional cross section in the Pacific at 160°W. There is also a significant contribution from region ⍀ 10 (not shown) in the top 100 m near the equator, with local maximums on both sides of the equator. 554 global maps of both ⌫ and ⌬ of the age distribution in order to characterize the time-scales associated with the ventilation of new water into the ocean interior.
We computed ⌫ and ⌬ directly without time stepping the transport equation by solving the stationary problems given by Eqs. (12) and (13) in section 2c. The solution was obtained by directly factoring the discretized matrix transport operator using a sparsematrix multifrontal approach (Amestroy et al. 2001) and then recursively back-solving to obtain the moments of the age distribution.
In Fig. 7, we illustrate spatial variations of the age distribution through plots of the zonally averaged ⌫(x) and ⌬(x). The upper-left panel shows that the oldest water (⌫ ϳ 900 yr) is situated in the northernmost part of the basin below 3000 m. The youngest waters are near the surface in the high-latitude regions close to deep convection regions and in the subtropical gyres where there is downward Ekman pumping. The standard deviation (upper-right panel) has a spatial pattern that is roughly similar to the mean age pattern. The broadest distributions (⌬ ϳ 650 yr) are in the northern part of the basin in the same region as the oldest water. The broad correlation between the (⌫) and (⌬) are a general feature of all the basins. The correlation is due to the fact that the width of the distribution is the result of the integrated effects of mixing since the water parcel left the surface. Since older water has had more time to experience mixing effects, it generally has a broader age distribution.
In the Pacific (middle panel), the oldest water with ⌫ ϳ 1600 yr is situated in the northern most part of the basin in the depth range between 2000 and 3000 m. The broadest distributions with ⌬ ϳ 1000 yr lie above the oldest water in the depth range between 1000 and 1500 m. In the Indian Ocean (bottom panel) the oldest water is again in the northern part of the basin below 1000 m where ⌫ exceeds 750 yr. The broadest age distributions in the Indian Ocean are located in a subsurface core centered at a depth of 750 m near 20°S where ⌬ reaches a maximum of 750 yr. Figure 8 shows horizontal maps of the mean lastpassage time (⌫) on the 792.5-, 1655-, and 3785-m depth levels. The top panel of Fig. 8 contours the mean on the 792.5-m depth level. Again we see how the Pacific water is older than the Atlantic water. The oldest water on the 792.5-m depth level has an age of approximately 900 yr and is situated in the western Pacific in a broad region directly north of the equator. In the Atlantic, the largest mean age is 300 yr and is situated in the eastern part of the basin near the equator. The middle panel shows the mean last-passage time on the 1655-m level. In the Atlantic there is a broad band of relatively young water along the western boundary that ages progressively as one moves southward. This age pattern is the result of the southward flowing deep western boundary current. In the South Pacific the mean last-passage time distribution is fairly zonal on the 1655-m level with the age increasing northward from 450 yr in the Southern Ocean to approximately 1250 yr near the equator. The oldest water on the 1655-m depth level is situated in the Northwest Pacific and in the region directly north of Indonesia with a mean age of more than 1400 yr. The bottom panel shows the mean first-passage time distribution on the 3785-m depth level. At this level, the mean age increases northward in both the Atlantic and Pacific reaching more than 950 yr in the North Atlantic and 1500 in the northeast Pacific.
Comparing the Atlantic and Pacific sections, we see that the deep water of the Pacific Ocean is much older than that of the Atlantic Ocean, consistent with the fact that there is no deep-water formation in the North Pacific. The surface waters near the equator are also much older in the Pacific Ocean. In the tropical Atlantic Ocean, the 200-yr ⌫ contour does not extend above 650 m, while in the tropical Pacific, it reaches up to a depth level of 150 m. The deep waters in the Indian Ocean are younger than those of the Pacific and generally older than those of the Atlantic.
The plots of ⌫ can also be compared with the idealage plots in England (1995). Some of the important differences between the two models include the higher vertical resolution, the different eddy mixing parameterization, as well as differences in the surface forcing. Our model has 29 levels in the vertical direction, whereas England's model had 12 levels. Our model uses GM and KPP mixing schemes, whereas England's model used constant horizontal and vertical eddy diffusion. Also the surface salinity in our model was forced by freshwater fluxes derived from an atmospheric climate model and was restored to the climatology with a 180-day time scale. In England's model salinity was restored to climatological values with a 50day time scale.
The major difference in the simulated age fields is that our model produced much older water in the deep North Atlantic than England's model; our model produced a zonally averaged mean age maximum of ϳ950 yr in the North Atlantic whereas England's model produced a maximum of ϳ500 yr. The age of the deep water in the Pacific was also older in our model, but the relative differences were not as big as in the Atlantic; our model produced a maximum zonally averaged mean age of ϳ1600 yr as compared with 1400 yr in England's model.
To determine which model result is more realistic we can use observed 14 C/C deficiencies in radiocarbon due to radioactive decay to estimate the time elapsed since a water parcel was last at the sea surface. For example, England (1995) compares his results with the 14 C/C sections measured as part of the Geochemical Ocean Sections Study (GEOSECS; e.g., Fig. 8 in England 1995).
One difficulty with such a comparison is that deepwater sources from the Northern and Southern Hemisphere have significantly different 14 C/C ratios. Broecker et al. (1998) point out that the southern end member has ⌬ 14 C ϭ Ϫ170% and the northern end member has ⌬ 14 C ϭ Ϫ67%. Consequently, the northto-south decrease in 14 C/C seen in the Atlantic Ocean is mainly the result of mixing between southern Circumpolar Water and North Atlantic Deep Water.
One way to handle this difficulty is to partition the interior water based on where it made its last contact with the surface so as to isolate the 14 C/C deficiencies due to radioactive decay. For the real ocean, Broecker et al. (1998) accomplished this by taking advantage of differences in the surface concentration of the quasiconservative tracer PO* 4 between the southern and northern water. We can use a similar idea to reconstruct the ⌬ 14 C distribution where, f ⍀ 1 and f ⍀ 5 are the fractions of northern and southern component waters and ⌫ ⍀ 1 and ⌫ ⍀ 5 are the means of the marginal age distributions for water that originated in regions ⍀ 1 and ⍀ 5 . See the appendix for a formulation of the mathematical equations used to compute the marginal age distributions, ⌫ ⍀ i . Figure 9 shows the model ⌬ 14 C distribution along the Atlantic GEOSECS section. Also shown in Fig. 9 is the fraction of northern (f ⍀ 1 ) and southern (f ⍀ 5 ) component waters along the section, together with the age of each component. The figure shows that the modeled carbon-14 deficiencies in the deep North Atlantic for our model are more than 2 times the measured value. [A plot of the Atlantic GEOSECS 14 C/C distribution can be found in England (1995) or in Broecker and Peng (1982).] The major cause for the discrepency is that Antarctic Bottom Water extends too far north. Estimates based on PO* 4 (Broecker et al. 1991) suggest that no less than 80% by volume of the water below 2000 m originated in the North Atlantic, but in our model 50%-75% of the water originated at the surface of the Southern Ocean. The larger initial 14 C/C deficiency of the southern water mass together with the longer transit time from the Southern Ocean to the North Atlantic both contribute to producing the excessively large 14 C/C deficiencies in our model. Exactly which feature of the model is the ultimate cause of this problem is beyond the scope of the present work-our purpose here is to point out the FIG. 9. Western Atlantic GEOSECS sections of (top left) f ⍀5 , (top right) ⌫ ⍀5 , (middle left) f ⍀1 , (middle right) ⌫ ⍀1 , (bottom left) ⌬ 14 C, and (bottom right) ⌫. The thick solid line is the f ⍀1 ϩ f ⍀5 ϭ 0.9 contour; above this contour, less than 90% of the water was formed in regions ⍀ 1 and ⍀ 2 . utility of using the last-passage time distribution for analyzing transport and ventilation in models. Deleersnijder et al. (2001) and Waugh et al. (2003) point out that the age derived from a decaying radioactive tracer is different from the mean age. For 14 C with a long half life when compared with the width of the age distribution the differences are small. Waugh et al. (2003) give the following approximation for the radio age in terms of the mean and standard deviation where Ϫ1 is the mean lifetime of a 14 C atom. For the North Atlantic where the maximum width is ⌬ ϭ 650 yr the resulting decrease in the radio age is approximately 26 yr for our model. In the North Pacific with a maximum width of ⌬ ϭ 1050 yr the decrease is about 67 yr. In both cases the correction due to the finite width of the age distribution is relatively small when dealing with the prebomb 14 C distribution. However, if one were dealing with the bomb 14 C distribution, which has relatively rapid atmospheric fluctuation, taking into account the finite width of the age distribution would be crucial for interpreting the 14 C/C deficiencies. The fi- nite width of the age distribution is also crucial for interpreting CFC ages (e.g., Waugh et al. 2003).
In Fig. 10 the zonally average of the mean (⌫ ) and width (⌬ ) of the first-passage time distribution is plotted. The spatial pattern of the mean first-passage time distribution is very different from that of the lastpassage time distribution, comparing with Fig. 7. The mean first-passage time has little latitudinal variation in comparison to that of the mean last-passage time. The largest mean first-passage times are in the bottom and exceed 1650 yr in each basin. As one moves up the water column, the mean first-passage time decreases and by the approximate middepth level near 2000 m it has decreased to roughly 600 yr before decreasing to zero at the surface. Mean first-passage times of less than 50 yr do not generally extend to regions deeper than 450 m except at high-latitude regions close to the regions of deep-water formation where they reach the bottom. The broadest first-passage time distributions (standard deviation Ͼ100 yr) are situated in the southern part of the Atlantic basin near the bottom (topright panel). The first-passage time distributions in the deep Pacific Ocean are generally less broad than in the Atlantic. This is consistent with the fact that the deep waters of the North Atlantic spread around the globe before making their first contact with the surface so that the deep Atlantic waters explore many more pathways to the surface than do the deep waters of the Pacific Ocean. In the deep Indian Ocean the firstpassage time distributions have widths that are of intermediate value between those in the deep Pacific and Atlantic oceans.
We conclude this section by summarizing the mean and standard deviation of the last and first-passage time distributions in tabular form. Table 2 gives the wholeocean mean and standard deviation of the last-and first-passage time in years partitioned by region where the last surface contact occurred. A map of the differ-ent regions is shown in Fig. 2, and the appendix gives a formulation of the equations that are solved to obtain the tabulated quantities. Table 2 shows that the North Atlantic (⍀ 1 ) and the Southern Ocean (⍀ 5 ) contribute on average the oldest water. This is as expected since the North Atlantic and the Southern Ocean are regions of deep-water formation, and the water that ventilates in those regions will remain below the surface for long time. Regions ⍀ 1 and ⍀ 5 also contribute most of the variance in the age distribution. This is consistent with the older mean age of water from these regionsmixing has been acting on the older water for a longer time. The mean age for the whole surface, 667 yr, is given in the last row of the table. Recall that the fraction of the whole ocean volume ventilated from each region is given in Table 1.
By contrast, the mean first-passage time is more uniform over the different regions. The exception being the regions in the Atlantic Ocean that have relatively short mean first-passage times. This reflects the fact the Atlantic is an exporter of deep water. The water that upwells in the Atlantic Ocean is therefore more likely to be upper ocean water with a shorter residence time in the ocean. The mean and standard deviation of the first-passage time to the whole surface are equal to the mean and standard deviation of the last-passage time distribution. This must be so, since in steady state, the fluxes of new and old water in and out of the ocean must balance in each age class.

The shape of last-passage time distributions
Typical shapes of the last-passage time distribution function, L ( lp | x), are shown in the right panel of Fig.  11 for four different positions, xЈ, in the interior of the ocean. They correspond to the locations on the 2352.5-m level indicated by the labels A, B, C, and D in the left panel of Fig. 11. The contours in the left panel indicate the volume fraction (f ⍀ 1 ) of water on the 2352.5-m level that was last at the surface in the North Atlantic. Nearly all of the water on the 2352.5-m level originates from either region ⍀ 1 (North Atlantic) or region ⍀ 5 (Southern Ocean). In fact, for x on the 2352.5-m level surface, indicating that nearly all of the water on this level originates from those two regions. The fluid volume in each of the grid boxes labeled A-D is composed of a mixture of fluid elements that were last at the surface at different times in the past. The plots in the right panel give the distribution of times when fluid elements made their last surface contact.
All four distributions have long exponentially decaying tails, with a universal e-folding decay time, 1 ϭ 797 TABLE 2. First and second column: mean (⌫ ⍀ i ) and standard deviation (⌬ ⍀ i ) of the volume-integrated last-passage time distribution from each ⍀ i region shown in Fig. 2. Third and fourth column: mean (⌫ ⍀ i ) and standard deviation (⌬ ⍀ i ) of the volumeintegrated first-passage time distribution to each of the ⍀ i region. The last row is for the first-and last-passage time to the whole surface ⍀ s . yr, given by the reciprocal of the lowest eigenvalue 2 of the transport operator. The universal exponential decay rate can be understood as follows: because the ocean has a finite volume, the transit time distribution can be written as a superposition of the transport operator's eigenmodes. Since all eigenmodes must decay, the contribution to the distribution from all but the lowest mode must eventually decay to expose only the most slowly decaying eigenmode. Close to the high-latitude deep-water formation regions the last-passage time distributions are peaked at young ages (locations A and D). As one moves away from the high-latitude regions of deep-water formation the last-passage time distributions become broader and peaked at progressively older ages (location B). In locations where there is confluence of different water masses, the last-passage time distributions can become multimodal (location C). At location C the distribution has two peaks corresponding to the different amount of time it takes for water from different sources to be transported there. As can be seen by the contours in the right panel of Fig. 11, approximately 35% of the water at location C originated from region ⍀ 1 . The other major contributor to the water at location C is water formed in the Southern Ocean, region ⍀ 5 . Since the time for fluid elements to be transported from North Atlantic Deep Water formation regions is much longer than the time for fluid elements to be transported from the Southern Ocean deep-water formation regions, the water formed in the North Atlantic is responsible for the second mode near 100 yr in the distribution. The water from the Weddell Sea is responsible for the first mode near 20 yr. Also, since the waters responsible for the older peak have had more time to mix with water of different ages, the older peaks are generally broader than the younger ones. Bimodal distributions similar to the one at location C exist in the regions close to the 0.3 and 0.4 contours in the left panel of Fig. 11 and also at similar locations on the levels above and below the 2352.5-m level in the South Atlantic and western Indian Ocean.

Discussion
The distribution of last-passage times from the surface provide a complete representation of the inte-2 The lowest eigenmode was computed using the inverse power method. grated transport to and from the surface. When this distribution is partitioned on where the last surface contact occurred, it provides a useful diagnostic tool for tracking the spreading of water masses from the surface formation location to any interior point.
An alternative representation of integrated transport is given by the distribution of first-passage times to the surface. When that distribution is partitioned based on where the first surface contact will occur, it gives information about when and where the current state of the ocean will become known to the atmosphere.
The problem of increasing atmospheric CO 2 concentrations provide nice applications for both first-and last-passage time distributions. Using last-passage time distributions, it is possible to inventory water masses into different age classes using atmospheric CO 2 concentrations from the time the water last visited the surface mixed layer. This can lead to better understanding of the anthropogenic CO 2 invasion into the deep ocean. For example, Hall et al. (2002) suggest fitting observed concentrations of several transient tracers to a parametric form of the last-passage time distribution and then to use the resulting distribution to propagate differences in dissolved inorganic carbon from preindustrial levels as a way of estimating anthropogenic carbon in the ocean. The form of the last-passage time distribution that they advocate using is an "inverse Gaussian" distribution. The advantage of this functional form is that it can be completely characterized by its first two moments, making it possible to use only two transient tracers to apply the method. The inverse Gaussian distribution may, however, provide a poor fit to the true distribution in some parts of the ocean. For example, our model results show that multimodal distributions are possible even for the case of steady transport. Thomas et al. (2001) use an approach similar to , but instead of using observed tracer data to determine the last-passage time distribution they use an idealized tracer in a global OGCM to compute the lastpassage time distribution (what they call age spectra).
The first-passage time distributions can be used to understand when the current state of the ocean will become known to the atmosphere. It has been suggested that CO 2 might be deliberately injected directly into the deep ocean as a way of mitigating atmospheric increases. It is therefore important to determine how effective such a strategy might be under different emission scenarios (Caldeira et al. 2001). The first-passage time distribution can determine when fluid elements in the deep ocean will be transported by ocean currents and mixing processes back to the surface mixed layer where their CO 2 will be able to leak back into the atmosphere.
The flux of the mean first-passage time through the surface in the adjoint flow is also a useful diagnostic since it tells us how much of the total ocean volume is ventilated through each surface region. Such a diagnos-tic is a fundamental property of the ocean circulation and is of direct interest to oceanographers. It also has useful applications for understanding the efficiency of the biological pump for sequestering CO 2 in the deep ocean.
The flux of the mean last-passage time (mean-age) through the surface tells us how much of the total ocean volume will become exposed to the atmosphere through each surface patch. This diagnostic can be useful for understanding the ultimate fate of different water masses.
In this paper, we have shown how to exploit the stationarity of the transport operator to compute the moments of the first-and last-passage time distributions by directly inverting the transport operator instead of time stepping the transport equation. The availability of fast matrix solvers for serial and distributed memory computers makes this an efficient way of obtaining integrated transport information directly from the differential form of the transport operator. In light of the fact that last-passage time distributions are increasingly being used in transient tracer studies, the matrix-inversion approach for computing the moments of the first-and last-passage time distribution might provide a useful tool for comparing the transport characteristics of different OGCM simulations. Further research is necessary however to determine the errors incurred by timeaveraging the transport operator.
The degree to which the transport operator can be approximated by a stationary transport operator will of course be problem dependent. In this paper we have used the time-averaged transport operator of a global ocean general circulation model with a coarse horizontal resolution. In this model of the ocean, the eddies are parameterized by eddy diffusion so that there is not any eddy variability. Time averaging was used to eliminate the seasonal cycle. We expect that the effects of the seasonal cycle should be confined mostly to the upper ocean. We plan to address the effects of the seasonal cycle on first-and last-passage time distributions in a follow up paper. It would also be interesting to study the effects of time dependence due to low-frequency variability in the strength of the circulation and in the end-of-winter mixed layer depth and in the intensity and location of deep-water formation. The first-and last-passage time diagnostic will also provide a useful diagnostic with which to compare parameterized eddydiffusion transport to ensembles of eddy-resolving simulations.
acknowledged. Useful comments from Juno Hsu and David Newman about the presentation and organization of the manuscript are also acknowledged.

Partitioned First-and Last-Passage Time Distribution
The mean and standard deviation of the marginal last-passage time distribution for the water fraction that made its last contact with the surface in ⍀ i will be denoted, respectively, by ⌫ ⍀ i and ⌬ ⍀ i ; that is, The mean and standard deviation of the marginal firstpassage time distribution for the water fraction that will make its first contact with the surface in ⍀ i will be denoted, respectively, by ⌫ ⍀ i and ⌬ ⍀ i :

͑A4͒
where ⌫ ⍀ i and ⌬ ⍀ i can be computed directly by solving stationary problems. To obtain the mean we solve subject to the surface boundary condititon ␥ ⍀ i ϭ 0, then set ⌫ ⍀ i ϭ ␥ ⍀ i /f ⍀ i . To obtain the standard deviation we solve subject to ␤ ⍀ i ϭ 0 at the surface, then set ⌬ 2 ⍀ i ϭ ␤ ⍀ i /f ⍀ i Ϫ ⌫ 2 ⍀ i . Here ⌫ ⍀ i and ⌬ ⍀ i can be obtained by solving equations similar to (A5) and (A6) in which the adjoint transport operator T † is used instead of T.