# Characterizing the chaotic nature of ocean ventilation

Graeme A. Macgilchrist, David P. Marshall, Helen L. Johnson, Camille Lique, Matthew Thomas

► To cite this version:

Graeme A. Macgilchrist, David P. Marshall, Helen L. Johnson, Camille Lique, Matthew Thomas. Characterizing the chaotic nature of ocean ventilation. Journal of Geophysical Research. Oceans, 2017, 122, pp.7577-7594. 10.1002/2017JC012875. insu-03682737

HAL Id: insu-03682737

https://insu.hal.science/insu-03682737

Submitted on 31 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

10.1002/2017JC012875

**Key Points:**

* A nondimensional "filamentation number" characterizes the chaotic nature of ventilation pathways

* Its magnitude in the subtropical North Atlantic implies a central role for mesoscale eddies

* The results have important implications for understanding the spreading of tracers in the ocean

**Correspondence to:**

G. MacGilchrist,
graeme.macgilchrist@univ.ox.ac.uk

**Citation:**

MacGilchrist, G. A., D. P. Marshall, H. L. Johnson, C. Lique, and M. Thomas (2017), Characterizing the chaotic nature of ocean ventilation, *J. Geophys. Res. Oceans*, 122, 7577–7594, doi:10.1002/2017JC012875.

Received 8 MAR 2017

Accepted 8 AUG 2017

Accepted article online 11 AUG 2017

Published online 19 SEP 2017

Characterizing the chaotic nature of ocean ventilation

Graeme A. MacGilchrist<fnref n="1" /> ID, David P. Marshall² ID, Helen L. Johnson<fnref n="1" /> ID, Camille Lique<fnref n="3" /> ID, and Matthew Thomas⁴

**Abstract** Ventilation of the upper ocean plays an important role in climate variability on interannual to decadal timescales by influencing the exchange of heat and carbon dioxide between the atmosphere and ocean. The turbulent nature of ocean circulation, manifest in a vigorous mesoscale eddy field, means that pathways of ventilation, once thought to be quasi-laminar, are in fact highly chaotic. We characterize the chaotic nature of ventilation pathways according to a nondimensional "filamentation number," which estimates the reduction in filament width of a ventilated fluid parcel due to mesoscale strain. In the subtropical North Atlantic of an eddy-permitting ocean model, the filamentation number is large everywhere across three upper ocean density surfaces—implying highly chaotic ventilation pathways—and increases with depth. By mapping surface ocean properties onto these density surfaces, we directly resolve the highly filamented structure and confirm that the filamentation number captures its spatial variability. These results have implications for the spreading of atmospherically-derived tracers into the ocean interior.

**Plain Language Summary** When water leaves the surface ocean and spreads into the ocean interior, it carries with it climatically important properties that have been exchanged with the overlying atmosphere, such as heat and carbon dioxide. It is likely that a significant part of this "ventilation" process is achieved by relatively small-scale (around 50–100 km) eddying motions, which are ubiquitous in the turbulent ocean, but this remains poorly understood and difficult to quantify. By drawing an analogy with the making of puff pastry - in which the baker thins the layers of dough by repeated stretching and folding - we propose a novel way of quantifying the role of eddying motions in ventilation. We evaluate the extent to which the eddying motions (the baker) generate thin filaments in a fluid parcel (the dough) in the ocean interior. This, in turn, indicates whether pathways of water from the surface ocean into the ocean interior are straightforward or "chaotic." In a numerical ocean simulation, we show that the latter is true - pathways are highly chaotic - supporting the case that eddying motions play an important role in the ventilation process.

1. Introduction

The process of ventilation sets the properties of the ocean subsurface through the subduction and spreading of mixed layer water. In the subtropical gyres it impacts thermocline stratification [Luyten et al., 1983] and the cycling of biogeochemical tracers [Karstensen et al., 2008] as well as affecting inter-annual and decadal climate variability through the exchange of heat [Deser et al., 1996; Williams et al., 2014] and carbon dioxide [Follows et al., 1996, 2002; Sabine et al., 2004] with the atmosphere.

Ocean circulation is highly turbulent, with the kinetic energy dominated by mesoscale eddies [Ferrari and Wunsch, 2009]. In the presence of strong eddy activity, the pathways of water parcels can be highly chaotic [Wiggins, 2005]. The chaotic nature of these pathways impacts the transport and distribution of tracers and the spatial scales over which regions of the ocean retain memory of their upstream properties.

Several studies indicate that the turbulent nature of ocean circulation is of central importance to thermocline ventilation [Robbins et al., 2000; Mecking et al., 2004; Booth and Kamenkovich, 2008; Koch-Larrouy et al., 2010; Abernathey and Ferreira, 2015; Kamenkovich et al., 2017]. Quantifying its impact, however, remains challenging due to the difficulty in separating out the effects of different spatial and temporal scales [Booth and Kamenkovich, 2008]. Consequently, analyses are often largely qualitative by necessity [e.g., Mecking

<footnote n="1">¹Department of Earth Sciences, University of Oxford, Oxford, UK, ²Department of Physics, University of Oxford, Oxford, UK,</footnote>

<footnote n="3">³Ifremer, Laboratoire d'Océanographie Physique et Spatiale, UMR6523, Brest, France, ⁴Department of Geology and Geophysics, Yale University, New Haven, Connecticut, USA</footnote>
et al., 2004; Kamenkovich et al., 2017). Notable exceptions to this are attempts to quantify, through a variety of means, a nondimensional Peclet number, which characterizes the relative importance of advection and mixing for tracer transport [e.g., Rhines and Schopp, 1991; Waugh et al., 2004; Shao et al., 2016].

In this contribution, we propose a novel, kinematic approach to characterize the chaotic nature of ocean ventilation pathways, borrowing concepts from the study of dynamical systems [Strogatz, 1994]. We define a nondimensional number that quantifies the reduction in the width of a fluid parcel since leaving the mixed layer. This indicates the extent to which trajectories from the mixed layer into the ocean interior are laminar or chaotic. The approach is complementary to tracer-derived estimates of ventilation pathways [e.g., Haine and Hall, 2002; Waugh et al., 2004; Primeau, 2005; Primeau and Holzer, 2006; Hall et al., 2007; Holzer et al., 2010], and explores the physical processes operating at their heart. Our analysis is focussed on the subtropical thermocline, a region of historical debate on the impact of ventilation and the time-varying flow on the large-scale circulation and stratification [Rhines and Young, 1982; Luyten et al., 1983].

Dynamical systems approaches such as that adopted in this study are often applied to questions of tracer stirring and transport in turbulent (and nonturbulent) geophysical fluids. These range from the investigation of Lagrangian flow maps that are mathematically similar to the advection-diffusion equation [e.g., Pierrehumbert, 1994; Antonsen et al., 1996; Neufeld et al., 2000; Tzella and Haynes, 2010] to evaluation of the dynamical system characteristics of observed or simulated currents and tracer fields [e.g., Abraham and Bowen, 2002; Waugh and Abraham, 2008; Haller, 2015; Mukhibi et al., 2016]. Oceanographic applications have tended to focus on the surface ocean (for a review, see Wiggins [2005]) and, to our knowledge, this study is the first to apply such an approach directly to the problem of ocean ventilation.

## 2. Theory

As an analogy to the effect of turbulent circulation on ocean ventilation, consider the making of puff pastry (Figure 1a). As the dough is repeatedly stretched and folded, the layers become increasingly thin. The change of the layer width, $\Delta x$, through time can be expressed as

$$ \frac{d\Delta x}{dt} = -\gamma\Delta x, \qquad (1) $$

where $\gamma$ is the strain rate—the rate at which the baker rolls the dough. Integrating (1), we find that the thinning of the layers occurs exponentially in time:

<figure><img src="image_3.png" /><figcaption>(b) Filamentation of fluid parcel</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Thinning of (a) layers of puff pastry (adapted from Strogatz [1994]) and (b) a patch of fluid at a circulation stagnation point, progressing in time from left to right. In both plots, the filament width is labeled $\Delta x$. In (a), the folding length is labeled $L$.</figcaption></figure>
$$ \Delta x(t) = \Delta x(0) e^{-\int_{0}^{t} \gamma dt} = \Delta x(0) e^{-\bar{\gamma}t}, \qquad (2) $$

where $\bar{\gamma}t$ is the time-mean strain rate. Thus the reduction in the width of the pastry layers is set by the time that the baker has been working for (t) and the vigour with which they roll the pastry ($\gamma$), which is averaged over time [Strogatz, 1994]. Importantly, it is the strain rate, $\gamma$, that impacts the layer width, while the horizontal scale of the pastry (denoted L in Figure 1a) is determined by the length scale of the folding.

Analogous to the stretching and folding of puff pastry is the horizontal drawing out and thinning of fluid parcels—so-called filamentation (Figure 1b). An infinitesimal fluid parcel of width $\Delta x$, aligned perpendicular to the axis of principal strain, will thin exponentially in time according to equation (2). In this case, the role of the baker is played by the circulation. For a two-dimensional, nondivergent flow, the strain rate, $\gamma$, is a function of the horizontal velocity gradients, given by

$$ 4\gamma^2 = \left( \frac{\partial u}{\partial x} - \frac{\partial v}{\partial y} \right)^2 + \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right)^2, \qquad (3) $$

where symbols have their usual meaning. As is the case for puff pastry, the filament width is affected by the strain rate, while the bounding horizontal scale of the fluid parcel is controlled by the characteristic length scale of the eddy field on which the fluid parcel is folded. Maxima in the strain rate are found at stagnation points in the flow. In a turbulent flow, where stagnation points are a common feature, the width of filaments will collapse to arbitrarily small scale in finite time.

Now consider a fluid parcel ventilated into the ocean subsurface. The time that the baker has been working for corresponds to the time since the fluid parcel left the ocean mixed layer—the ventilation timescale. The vigour with which the baker rolls the pastry is the strain rate of the circulation, which varies in space and time along the parcel’s Lagrangian trajectory. Then, from (2), the reduction in filament width is a function of the ventilation timescale ($\tau_{vent}$) and the mean strain rate along the Lagrangian trajectory ($\bar{\gamma}^L$; we use the superscript L to distinguish this Lagrangian mean from the Eulerian time-mean strain rate of the puff pastry evaluated at a stationary point). Expressing the latter as a Lagrangian strain timescale ($\bar{\tau}_{strain}^L = 1/\bar{\gamma}^L$), we can define a nondimensional “filamentation number” as the ratio of these two timescales:

$$ F = \frac{\tau_{vent}}{\bar{\tau}_{strain}^L}. \qquad (4) $$

The filament width of a fluid parcel will have been reduced by a factor of order $e^{-F}$ since leaving the mixed layer (equation (2)). For example, in a region where $F = 4$, a ventilated fluid patch with an initial horizontal scale of $\Delta x(0) = 100$ km will have formed filaments of roughly 2 km width (notwithstanding spatial scale considerations outlined below).

In other geophysical fluid applications of concepts from dynamical systems, the filamentation of fluid parcels is characterized by the “Lyapunov exponents” of the system. The finite time Lyapunov exponent (FTLE) is the time mean strain rate following a trajectory. Thus, for individual trajectories, the FTLE is identical to the Lagrangian strain rate defined above. Finite time Lyapunov exponents are evaluated from trajectories initialized across a domain and evolved forward or backward for a fixed length of time [e.g., Abraham and Bowen, 2002]. Their evolution and distribution as the fixed time increases, and their long time behaviour, has been shown to provide information on the stirring of passive tracers within the domain [Antonsen et al., 1996; Abraham and Bowen, 2002]. In the present work, the finite time of integration for each trajectory is the time since the parcel was last in the surface mixed layer, i.e., the ventilation timescale. As such, it is unique for each trajectory and each location within the domain. It is for this reason that we identify the Lagrangian strain rate as a distinct measure from the more ubiquitous FTLE.

The formulation of the nondimensional number F assumes that the width of the fluid parcel is infinitesimally small. For F to characterize the filamentation of a finite-size fluid patch, the width of the patch must be less than the dominant kinetic energy-containing spatial scale of the circulation. That is to say the fluid
<figure><img src="image_5.png" /></figure>

patch must be within a “nonlocal” dynamical regime, in which the kinetic energy spectrum of the circula-
tion, and therefore the strain timescale, is dominated by the larger scales [Batchelor, 1959; Bennett, 1984;
Keating et al., 2011]. To understand this, consider a circulation regime for which the dynamics are nonlo-
cal below 100 km. A ventilated fluid patch with initial horizontal scale of 50 km will be persistently
thinned on a strain timescale associated with the nonlocal regime. This timescale is determined by the
circulation at the largest scales of the nonlocal regime—100 km. In this case, *F* will correctly characterize
the reduction in filament width. On the other hand, a ventilated fluid patch with an initial horizontal
scale of 500 km—within a “local” regime—will experience a different strain (i.e., that associated with a
separation of 500 km), the timescale of which is likely to be longer than that within the nonlocal range.
Consequently, *F* will generally overestimate the extent of filamentation of such a fluid patch. Returning
briefly to the puff pastry analogy, the baker cannot efficiently roll the dough if its size exceeds their
own.

The puff pastry analogy is commonly employed to explain the mechanism by which dynamical systems
develop sensitivity to initial conditions, through the stretching and folding of phase space [Strogatz, 1994].
Imagine a small drop of food colouring in the pastry dough (Figure 1a), representing a cluster of nearby ini-
tial conditions in phase space. After a number of iterations of stretching and folding, the dye will have
spread throughout the pastry, representing the range of possible trajectories of the initially close initial con-
ditions. As the layers thin, trajectories in the system become increasingly sensitive to their initial position —
the system becomes increasingly chaotic.

By the same analogy, the value of *F* in the ocean subsurface characterizes the chaotic nature of ventilation
pathways. This is illustrated schematically in Figure 2, where we consider the ventilation pathways that
have been followed by adjacent particles on a density surface. If a region is characterized by a small value
of *F*, implying that it is ventilated on a short timescale (τ<sub>vent</sub> small) or the flow is predominantly laminar
(τ<sub>strain</sub><sup>L</sup> large), adjacent particles will have left the mixed layer at similar times and locations (blue path-
ways). In the context of a passive tracer, its distribution on leaving the mixed layer will be conserved in
the ocean interior. If, on the other hand, the region is characterized by large *F*, suggesting it is ventilated
on a long timescale (τ<sub>vent</sub> large) or the flow is highly turbulent (τ<sub>strain</sub><sup>L</sup> small), the particles are more likely
to have followed divergent pathways and have left the mixed layer at a range of locations and times (red
pathways). In this case, the distribution of a passive tracer will be highly filamented, reflecting these diver-
gent pathways.

In the following, we explore the concept of the filamentation number, evaluating its magnitude within the
subtropical thermocline of a numerical ocean circulation model and comparing it to the ventilation path-
ways in this region.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Schematic illustration of laminar (blue) and chaotic (red) ventilation pathways of adjacent particles on a density surface. The particles are adjacent at time *t* and the *i*-th particle left the mixed layer at time *t*—*Δt*<sub>*i*</sub>.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

## 3. Methods

### 3.1. Model Details and Lagrangian Trajectory Calculation

Our analysis utilizes the 5-day mean output from a hindcast simulation performed with the eddy-permitting ORCA025 coupled ocean/sea-ice model configuration developed by the European Drakkar project [Barnier et al., 2006], which is an implementation of the NEMO model [Madec, 2008]. The configuration uses an isotropic tripolar grid covering the global ocean with a resolution of 1/4° (~27.75 km at the equator, refined at higher latitudes) and 75 uneven vertical levels. The simulation is interannual and runs from 1958 to 2010 with no spin-up. The forcing data set is the Drakkar forcing set 5, which is an updated version of the fields described in Brodeau et al. [2010]. Further details and observational validation can be found in Grégorio et al. [2015] and Barrier et al. [2015].

We also make use of the particle tracking tool, ARIANE (http://www.univ-brest.fr/lpo/ariane/) to calculate Lagrangian trajectories from the hindcast velocity fields. A description of the algorithm is given by Döös [1995] and Blanke and Raynaud [1997]. The algorithm makes use of the fact that, for a numerical simulation constructed on a C-grid, the stream function between model grid faces can be determined analytically. This relies only on the assumption that velocities vary linearly between grid faces. Since we compute our Lagrangian trajectories from 5 day mean fields, we neglect velocity and tracer fluctuations, and their correlations, on timescales shorter than 5 days.

### 3.2. Timescale Calculations

Both the ventilation and strain timescales are evaluated using Lagrangian trajectories. Particles are initialized across potential density surfaces at a resolution of 25 per grid cell (~5 km apart) in March 2010. They are tracked backward in time until they intercept the mixed layer (which we define by a density within 0.01 kg m⁻³ of the surface density, following Thomas et al. [2015]). Along the trajectory, the local horizontal strain rate, determined from the Eulerian velocity field, is linearly interpolated onto the trajectory location, providing a time series of strain rate for each trajectory. We use the horizontal strain rate under the assumption that flow in the subsurface is predominantly along density surfaces. The error associated with the density surface slope is negligible. To determine $\bar{\tau}_{\text{strain}}^L$, we calculate the inverse of the mean strain rate along each trajectory, and grid this on the particle's initial location on the density surface. For $\tau_{\text{vent}}$, the age of each particle (the time between when it left the mixed layer and March 2010) is gridded on its initial location on the density surface. Median values of $\bar{\tau}_{\text{strain}}^L$ and $\tau_{\text{vent}}$ are determined within each model grid cell (25 particles), and the resulting fields are further smoothed using a two-dimensional, 24 grid cell, spatial moving median. The result is a smoothed map of the ventilation timescale and the Lagrangian strain timescale across the density surface.

### 3.3. Structure Function Calculation

The capacity for F to characterize filamentation of finite-size fluid patches depends on the details of the oceanic energy spectrum, specifically the spatial scales over which the dynamics are nonlocal (section 2). We determine this by analysis of the Eulerian velocity field, from which we evaluate the variance of separation velocities — the second order structure function — given by:

$$S_2(\mathbf{r}) = \langle [u(\mathbf{x}+\mathbf{r}) - u(\mathbf{x})]^2 \rangle, \quad (5)$$

where r is separation and angular brackets denote the time average. $S_2$ can be evaluated for velocities and separations in both zonal and meridional directions. From dimensional analysis, the dependence of $S_2$ on $|\mathbf{r}|$ indicates the dynamical regime of the velocity field [LaCasce, 2002]. Where $S_2 \propto |\mathbf{r}|^{2/3}$, the dynamics are local, equivalent to an energy spectral slope of -5/3. Where $S_2 \propto |\mathbf{r}|^2$, the dynamics are nonlocal, equivalent to an energy spectral slope of -3. In our analysis, second order structure functions are evaluated using horizontal velocities interpolated onto density surfaces in the thermocline.

## 4. Results

### 4.1. Quantifying the Filamentation Number

Our analysis considers three potential density surfaces in the North Atlantic subtropical thermocline (Figure 3). These surfaces, $\sigma_0=26 \text{ kg m}^{-3}$ (~100–300 m depth), $\sigma_0=26.5 \text{ kg m}^{-3}$ (~200–600 m depth) and $\sigma_0=27 \text{ kg m}^{-3}$ (~300–700 m depth), were chosen to sample a range of dynamically distinct regimes in
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Depth of three potential density surfaces in March 2010: (a) $\sigma_{\theta}=26 \text{ kgm}^{-3}$, (b) $\sigma_{\theta}=26.5 \text{ kgm}^{-3}$ and (c) $\sigma_{\theta}=27 \text{ kgm}^{-3}$.</figcaption></figure>

the subtropical thermocline [Samelson and Vallis, 1997]. The shallower two surfaces both outcrop in the downwelling region of the subtropical gyre, while the deeper surface largely outcrops to the north of the Gulf Stream and North Atlantic Current, except at its eastern extent. For the rest of this paper, we refer to these surfaces as $\sigma_{26}$, $\sigma_{26.5}$ and $\sigma_{27}$.

Ventilation timescales across these three surfaces are shown in Figure 4. Across most of the subtropical gyre, $\sigma_{26}$ is ventilated on timescales less than 5 years (Figure 4a). Recently ventilated water, with $\tau_{\text{vent}}$ less than 2 years, stretches across half of the basin. The intermediate density surface, $\sigma_{26.5}$, is mostly ventilated on timescales between 5 and 15 years (Figure 4b). The youngest ventilated fluid is restricted to a relatively small (by comparison with that on the shallow surface) region east of 20°W. On both of these surfaces, longer ventilation timescales are seen in the so-called "shadow zone"—a region of stagnant mean flow off the west coast of Africa that is not directly ventilated by the wind-driven gyre [Rhines and Young, 1982; Luyten

<figure><img src="image_10.png" /><figcaption>**Figure 4.** The ventilation timescale, $\tau_{\text{vent}}$, in years. Values are averaged within each model grid cell and further smoothed by averaging with the adjacent 24 grid cells. Locations in black were unventilated during the model run. Locations in light grey are where more than half of the particles in the region originate outside the experimental domain, thus no ventilation timescale information is obtained.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 5.** (a, b, c) The Eulerian strain timescale in years, evaluated at each horizontal location on the three density surfaces ($\tau_{\text{strain}}^E$), averaged over the last 30 years of the numerical simulation. (d,e,f) The Lagrangian strain timescale ($\tau_{\text{strain}}^L$) in years, evaluated as the inverse of the mean strain rate experienced along trajectories that are present on the density surfaces in March 2010, and that left the surface ocean within the time of the model simulation. Values are averaged within each model grid cell and further smoothed by averaging with the adjacent 24 grid cells. Locations in black were unventilated during the model run. Locations in light grey are where more than half of the particles in the region originate outside the experimental domain, thus the full strain history could not be determined. Note the logarithmic color scale in both Figures 5a and 5b. (g, h, i) The ratio of the Lagrangian and Eulerian strain timescales, presented on a log scale.</figcaption></figure>

et al., 1983). The deeper isopycnal, $\sigma_{27}$, is ventilated on considerably longer timescales (Figure 4c). A gradient exists between comparatively short timescales on the order of 1–10 years near the density surface outcrop and longer timescales of 50 years and more to the south west of the gyre. There is a clear signal of long ventilation timescales within the Gulf Stream and its extension, implying that water in this region is recirculated from the poorly ventilated southwest gyre interior. The distribution of $\tau_{\text{vent}}$ on these surfaces is broadly consistent with observations of tracer-derived ages [Jenkins, 1988] and previous model results [Williams et al., 1995].

Before presenting the Lagrangian strain timescale, $\tau_{\text{strain}}^L$, we consider the pattern of the time-averaged, local, Eulerian strain timescale, $\tau_{\text{strain}}^E$, across the density surfaces (Figures 5a–5c). This is determined at each grid cell location using the horizontal Eulerian velocity field interpolated onto the potential density surfaces
<figure><img src="image_13.png" /></figure>

and averaged over the last 30 years of the simulation. The pattern of $\tau_{\text{strain}}^E$ on each density surface provides an overall impression of where trajectories will experience high or low strain. To the best of our knowledge, this informative dynamical quantity has not previously been mapped in a realistic ocean circulation model. The magnitude of $\tau_{\text{strain}}^E$ is similar on the three density surfaces, increasing slightly with depth in the gyre interior and equatorial region. It ranges from days in energetic regions such as the Gulf Stream and equatorial band, to months in the gyre interior and shadow zone, and is everywhere less than a year. Within the comparatively quiescent gyre interior, there remains considerable spatial structure, with $\tau_{\text{strain}}^E$ increasing south of the Gulf Stream and decreasing again toward the gyre's southern extent. From observations of a tracer released at 25.6°N, 28.3°W, Ledwell et al. [1998] estimate a strain timescale of ~0.1 years at a potential density of 26.75 kgm⁻³, consistent with the time-average value in this region presented here. On $\sigma_{26.5}$ and $\sigma_{27}$, a band of shorter timescales on the east of the basin at around 35°N signals the impact of the Azores current and possibly its interaction with Mediterranean water eddies [Bashmachnikov et al., 2015]. There is a notable increase in $\tau_{\text{strain}}^E$ on all three density surfaces as the boundary between the gyre circulation and shadow zone is crossed.

The Lagrangian strain timescale ($\bar{\tau}_{\text{strain}}^L$; Figures 5d–5f) is the inverse of the mean strain rate experienced by the fluid parcels along their trajectories from the mixed layer to each of the three density surfaces. Thus, it integrates the local strain rate as trajectories move around the subsurface. Consequently, and consistent with finite time Lyapunov exponents [Antonsen et al., 1996; Abraham and Bowen, 2002; Waugh et al., 2006], values of $\tau_{\text{strain}}^L$ should converge at long times. The time for convergence of FTLEs in the surface ocean (~1 year) [Abraham and Bowen, 2002; Waugh et al., 2006] is considerably less than the ventilation timescale across most of the thermocline in our simulation (Figure 4). It is therefore unsurprising that the distribution of $\bar{\tau}_{\text{strain}}^L$ is relatively uniform compared to $\tau_{\text{strain}}^E$. The remaining basin-scale variation in $\bar{\tau}_{\text{strain}}^L$, and the variation between surfaces, must be due to the large-scale pathways of the trajectories. On $\sigma_{26}$, where trajectories have been comparatively recently ventilated, $\bar{\tau}_{\text{strain}}^L$ retains some of the spatial structure of the $\tau_{\text{strain}}^E$ distribution, with longer timescales in the central gyre. In the older waters present on deeper surfaces, trajectories are more likely to have passed through regions of short local strain timescale, particularly given the outcrop location of $\sigma_{26.5}$ and $\sigma_{27}$ to the north of the Azores Current and North Atlantic Current (Figures 5a–5c). This results in a slight decrease in $\bar{\tau}_{\text{strain}}^L$ with depth. By contrast, $\tau_{\text{strain}}^E$ shows a marginal overall increase with depth. On both of the deeper surfaces, $\bar{\tau}_{\text{strain}}^L$ increases toward the south of the gyre. A large number of particles at the southern extent originate in the ventilation region to the north east of the density surfaces. Their trajectories, therefore, avoid passage through the most turbulent region with shortest strain timescale, i.e., the Gulf Stream and North Atlantic Current. It is apparent, therefore, that the distribution of $\bar{\tau}_{\text{strain}}^L$ is strongly influenced by the overall pathways of ventilation, and the strain history along these pathways, rather than being entirely set by the local values of Eulerian strain on each surface.

An Eulerian strain timescale is more easily evaluated than a Lagrangian strain timescale (both in a numerical model and in the ocean). Consequently, it is instructive to understand the extent to which $\tau_{\text{strain}}^E$ differs from $\bar{\tau}_{\text{strain}}^L$, and whether the former may be adequate for the purpose of scaling analysis. To this end, their ratio on each density surface is presented in Figures 5g–5i. On each density surface, differences between $\bar{\tau}_{\text{strain}}^L$ and $\tau_{\text{strain}}^E$ exist in the quiescent regions (the gyre interior) and energetic regions (such as the Gulf Stream), where the Eulerian value is longer and shorter respectively. Trajectories in these regions have previously experienced longer/shorter strain timescales than the short/long local $\bar{\tau}_{\text{strain}}^L$, and these are incorporated into the calculation of $\tau_{\text{strain}}^E$. Nevertheless, the difference between the Eulerian and Lagrangian strain timescale is much less than an order of magnitude across almost all of the area of the density surfaces, implying that an Eulerian measure may be broadly adequate for the scaling analysis in this paper.

The filamentation number, the ratio of ventilation and Lagrangian strain timescales, is shown in Figure 6. It is large ($F \gg 1$) throughout the thermocline, except in a narrow strip adjacent to the winter outcrops, and increases markedly with depth. The magnitude of F is a consequence of the considerably longer timescale on which the ocean subsurface is ventilated compared to the timescale over which it is strained. Owing to the relative uniformity of $\bar{\tau}_{\text{strain}}^L$, the spatial distribution of F on each surface and with depth is almost entirely set by that of $\tau_{\text{vent}}$. This is highlighted by plotting F against $\tau_{\text{vent}}$ for all density surfaces in Figure 7. As such, the lowest values of F are seen at the eastern end of the outcrops, where waters are recently
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 6.** The filamentation number, $F = \tau_{vent} / \tau_{strain}^L$, on the three density surfaces. Locations in black were unventilated during the model run. Locations in light grey are where more than half of the particles in the region originate outside the experimental domain, thus neither the ventilation time or the strain history could be determined. Green, dotted lines show the location of high resolution mappings relevant to section 4.3.</figcaption></figure>

ventilated, and increase westward. On the shallower surface, $F$ increases from around 10 to more than 100 across the basin, while on the intermediate surface $F$ approaches 300 toward the Gulf Stream region. Across the deeper surface, $F$ exceeds 300 everywhere except in ventilation regions off the coast of Portugal and northward of the North Atlantic Current.

As outlined in section 2, $F$ characterizes filamentation of fluid patches with an initial width within the nonlocal dynamical regime. The spatial range of the nonlocal regime can be determined by calculating the variance of separation velocities of the Eulerian velocity field—so-called structure functions (section 3.3). In this numerical simulation, the structure function is proportional to the square of the separation on scales from the grid scale up to around 100 km, as illustrated for a range of locations across the three density surfaces in Figure 8. This implies that the dynamics are nonlocal within this range, and that $F$ characterizes filamentation for fluid patches with horizontal scales up to at least 100 km.

<figure><img src="image_16.png" /><figcaption>**Figure 7.** Two-dimensional histogram of ventilation timescale ($\tau_{vent}$) and filamentation number ($F$), where the data are from all three density surfaces. Note the logarithmic color scale.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 8.** Second-order structure functions ($S_2$) at a range of locations (denoted by the line colors) across each of the three density surfaces (denoted by the line styles). See section 3.3 for details of the calculation. A line corresponding to $S_2 \propto r^2$ is also shown.</figcaption></figure>

## 4.2. Mapping Chaotic Pathways

The magnitude of the filamentation number characterizes the chaotic nature of ventilation pathways. We now determine these pathways directly by evaluating Lagrangian trajectories backward in time from an array of points on the three density surfaces to the time and location where they intercept the surface mixed layer. This procedure defines a mapping from each initial position on the density surface to the year and location of ventilation, portrayed schematically in Figure 2. The year and location act as nondiffused passive tracers, their distribution set on subduction from the mixed layer and stirred in the subsurface by the circulation. Akin to the mapping of solutions in dynamical systems [see e.g., Strogatz, 1994; Pierrehumbert, 1994; Neufeld et al., 2000], the pattern reveals the chaotic nature of the ventilation pathways.

The mappings are produced in March 2010 across each density surface (Figures 9a–9c, 10a, and 10b). Additionally, we zoom in on interesting regions (Figures 9d–9f, 10c, and 10d). While the horizontal resolution of the velocity field is 1/4°, the mapping is calculated at higher spatial resolution, which is possible due to the interpolation of velocities between grid cells in Lagrangian analysis (section 3.1). The mapping is purely advective and isolates the role of the model’s resolved scales in producing filamentary structure. In the ocean, the filamentary structure produced by the larger-scale velocity field ceases when local dynamics, such as submesoscale motions, become dominant [Smith and Ferrari, 2009; Callies and Ferrari, 2013].

Most immediately striking in Figure 9 is the contrast between the mappings on different density surfaces (Figures 9a–9c), which becomes increasingly noisy and chaotic with depth. The proportion of the density surfaces occupied by coherent regions of constant ventilation year (a “homogenous” mapping) is sharply reduced between $\sigma_{26}$, where recently ventilated “2009” water occupies the eastern third of the basin, and $\sigma_{26.5}$, where it is restricted to the east of 20°W. On both of these surfaces, the mapping becomes increasingly filamented and chaotic moving westward away from these regions. On $\sigma_{27}$, similar homogenous regions are absent, and the mapping is comparatively chaotic across the entire surface. Such variation between surfaces is consistent with our expectations based on a filamentation number that increases with depth (Figure 6).

Considering initially $\sigma_{26}$, the homogenous nature of the mapping in the eastern half of the basin is perhaps surprising considering the values of F in this region (between 10 and 50) from which one would expect to find smaller scale filaments present. That this filamentation is not present is due to the large spatial scale occupied by recently ventilated 2009 water, which exceeds the scale over which filamentation characterized by F occurs (~100 km, section 4.1, Figure 8). Consequently, filamentation of 2009 water is confined to its outer edges. This can be seen in the centre of the gyre and, most strikingly, at the boundary between the
<figure><img src="image_19.png" /></figure>

Year of ventilation

<figure><img src="image_20.png" /><figcaption>**Figure 9.** Ventilation year of particles on (a, d) $\sigma_{\theta}=26 \text{ kgm}^{-3}$, (b, e) $\sigma_{\theta}=26.5 \text{ kgm}^{-3}$ and (c, f) $\sigma_{\theta}=27 \text{ kgm}^{-3}$ at the end of March 2010. (a, b, c) Entire density surface, 25 particles per model grid cell. (c, d, e) An approximately $10^{10} \times 10^{10}$ box centred on (c) 25°W 20°N, (d) 30°W 30°N and (e) 15.5°W 39.5°N, 625 particles per model grid cell, ~1 km resolution. Note that the color scale is nonlinear, designed to emphasize gradients in the most recent two decades. Black: particle unventilated during model run. Light grey: Particle originated outside model domain, thus no ventilation year obtained.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 10.** Longitudinal band from which particles were ventilated on (a, c) $\sigma_{\theta}=26 \text{ kgm}^{-3}$ and (b, d) $\sigma_{\theta}=26.5 \text{ kgm}^{-3}$ at the end of March 2010. (a, b) Entire density surface, 25 particles per model grid cell, and (c, d) approximately $10^{\circ} \times 10^{\circ}$ boxes centred on 38° W 22° N, 625 particles per model grid cell, ~1 km resolution. Black: particle unventilated during model run. Light grey: Particle originated outside model domain, thus no ventilation longitude obtained.</figcaption></figure>

gyre interior and shadow zone to the south (Figures 9a and 9d). The homogenous mapping of recently ventilated water contrasts starkly with the more filamented, chaotic shadow zone, where adjacent water parcels are often ventilated decades apart. Long, thin filaments of the shadow zone water are drawn out into the newly ventilated gyre interior, and vice versa (Figure 9d). This eddy-driven exchange is hypothesized to be crucial for the ventilation of shadow zones [Bettencourt et al., 2015].

Filamentation occuring within the recently ventilated 2009 water on $\sigma_{26}$ can be visualized by mapping ventilation longitude rather than ventilation year (Figure 10). Given the predominantly zonally-aligned nature
<figure><img src="image_23.png" /></figure>

of the density surface outcrop, and by construction of the color scale, the maximum initial scale of venti-
lated water in this mapping is 5° longitude (≈ 500 km), and therefore closer to the scales at which the fila-
mentation characterized by *F* occurs. Consequently, filamentation occurs within a short distance of the
outcrop, including within the recently ventilated 2009 water (Figures 10a and 10c). The mapping of ventila-
tion longitude is also indicative of the large horizontal scale over which stirring and filamentation mixes
fluid. Regions of the 10° by 10° box on σ₂₆ (Figure 10c) show trajectories that were ventilated from across a
45° longitude band, within a year of ventilation.

On σ₂₆.₅, recently ventilated water is confined to a comparatively small region at the eastern side of the basin (Figure 9b). The mapping becomes increasingly chaotic within a short distance to the west, where filaments of young and old water are stirred together (Figure 9e). The contrasting nature of the ventilation pathways on σ₂₆.₅ and the shallower surface, σ₂₆, is emphasized by considering the mapping of ventilation longitude (Figure 10). The high-resolution mappings in Figures 10c and 10d are positioned at the same horizontal location. Consequently, they exhibit the same structure of recent mesoscale stirring on both surfaces. On σ₂₆.₅ the stirring is taking place between already highly filamented water, due to the greater age of the water on this surface and the strain it has experienced since ventilation — characterized by its greater value of *F*.

The mapping of ventilation longitude gives an impression of the large-scale ventilation pathways on each of these shallower surfaces. On σ₂₆, despite filamentation taking place, there remains a strong signature of the gyre circulation (Figure 10a). There is some evidence for the impact of the large-scale gyre circulation on σ₂₆.₅, where fluid ventilated at the eastern edge of the outcrop appears to be preferentially advected around the southern edge of the gyre, though there is evidently significant filamentation and stirring along the way (Figure 10b). Both surfaces exhibit increasingly filamented structure toward the western “pool” region, in which eddy stirring is expected to be dominant [Rhines and Young, 1982]. On both surfaces, there appears to be eastward movement of water immediately south of the outcrop, which appears to be due to eddy stirring and, on σ₂₆.₅, advection by the Azores current. Finally, the mapping of ventilation year (Figure 9) gives the impression that, on both of these surfaces, water is ventilated in the east and ages as it moves westward. In contrast to this view, the western pool region on both density surfaces is at least partly ventilated locally in the west rather than being composed of waters of eastern origin (Figures 10a and 10b). This raises the possibility that the age distribution on these surfaces is, in part, due to stirring between a region in the east that is locally ventilated each year and a region in the west that is ventilated locally, but on much longer long-timescales.

Across the entirety of the deepest surface, σ₂₇, the mapping of ventilation year is highly chaotic, consistent with the large filamentation number, and suggesting that the ventilation pathways of this density surface are dominated by eddy-stirring (Figure 9c). There persists a large-scale gradient in ventilation year from the northeast to the southwest of the basin, often invoked to hypothesize a dominant role for the gyre circulation [Williams, 1991; Jenkins, 1988]. However, such a gradient could feasibly be established by eddy-stirring alone, from the ventilation region in the north east, without any role for the time-mean flow. Such a scenario would be facilitated by the anisotropic nature of eddy-stirring in the North Atlantic, which appears to broadly align with streamlines of the large-scale flow [Rypina et al., 2012].

Even close to the ventilation region in the north east, the mapping is remarkably chaotic, with particles ini-
tialized side-by-side often found to have been ventilated decades apart (Figure 9f). Recently ventilated
water, from 2009, is filamented and stirred into older water close to the density surface outcrop. Within one
year, the filaments of 2009 water have been reduced to extremely small scales (Figure 9f), emphasizing the
exponential rate at which this filamentation occurs. The short timescale of filamentation is similarly appar-
ent on the shallowest surface, σ₂₆, at the shadow zone boundary (Figure 9d) and in the gyre interior (Figure
10c). In both of these regions, the observed filaments have been generated in less than one year.

4.3. Validating the Filamentation Number

The theory presented in section 2 indicates that the filamentation number, *F*, characterizes the filament width (and therefore the chaotic nature) of the mappings in Figures 9 and 10. While obtaining a direct measure of filament width is not straightforward, we can investigate properties of the tracer distribution that may be expected to be related, such as power spectra and tracer gradients. Given the magnitude of *F*, we expect the filament width to be extremely small — orders of magnitude smaller than the model grid scale. To resolve as much of the filamentation as possible, therefore, we produce a mapping at very high spatial
<figure><img src="image_24.png" /></figure>

resolution (approximately 1 m) along zonal lines on each of the three surfaces. The location of these zonal lines is shown in Figure 6, from which it can be seen that they sample a range of values of F.

The power spectra of passive tracers in the nonlocal regime are predicted to exhibit power law behaviour, with the power varying as $k^\alpha$ where k is the spatial wavenumber [Batchelor, 1959]. Pierrehumbert [1994] found that a slope of $\alpha=-1$ (consistent with the prediction of Batchelor [1959]) emerges in the nondiffusive limit of a dynamical system map developed to simulate chaotic advection. They also note that a flattening of the spectrum occurs when the expected filament width of the tracer is much less than the resolution of the mapping as a result of aliasing from undersampling. That is to say, the presence of increasingly small scales impacts the spectral slope in the resolved range. It logically follows, and indeed Pierrehumbert [1994] found, that the smaller the expected filament scale, the flatter the spectral slope. While this aliasing obscures the signal that Pierrehumbert [1994] were investigating, in our case we would expect its effect to be directly reflected by the filamentation number — a larger value of F implies the presence of smaller scales and consequently a flatter spectral slope. The spatial power spectrum of the ventilation longitude is evaluated within each grid cell (2.5×10⁴ data points) using Welch's power spectral method [Welch, 1967]. Grid cells are grouped according to their value of F, and the average spatial power spectrum in different ranges of F is determined.

The results of the spectral analysis can be seen in Figure 11a. We only show averages for regions with F less than 500. For F less than 25, the spectral slope is slightly steeper than -1. This is likely to be due to the fact that the forcing of tracer variance — in this case ventilation — is more similar to a "resetting" forcing [Pierrehumbert, 1994; Neufeld et al., 2000] than an "additive" forcing, adopted in the theory of Batchelor [1959]. As we consider regions of larger F, the slope of the average spectrum become increasingly shallow, consistent with smaller filament width, and a more chaotic mapping in these regions.

The filament width should also be related to the gradient of the tracer in the mapping. The tracer gradient, $\nabla C$, scales as $\Delta C/\Delta x$, where $\Delta x$ is the filament width and $\Delta C$ is a typical tracer difference between adjacent filaments. So as $\Delta x$ decreases, $\nabla C$ increases. The PDFs of the squared gradient of the ventilation longitude are shown in Figure 11, where the x-axis is on a natural logarithm scale. As for the power spectra, PDFs are grouped and averaged within bands of F. For the ventilation longitude, the relationship between $\nabla C$ and F is complicated by the fact that $\Delta C$ will also vary depending on the locations where, and scales across which, subduction occurs, as well as the spatial scales of stirring in the subsurface. Close to a ventilation region, $\Delta C$ will tend to be small, with water parcels mostly ventilated locally. Deeper in the subsurface, where parcels have originated from a wider range of longitudes, its typical value will be larger. We correct for this by normalizing the tracer gradient by the standard deviation of the tracer concentration in each grid cell.

Figure 11b shows that at lower values of F, normalized tracer gradients tend to be smaller, consistent with our expectations. If the water were subducted uniformly across the density surface outcrop and not stirred

<figure><img src="image_25.png" /><figcaption>**Figure 11.** (a) Average spatial power spectrum and (b) PDF of gradients (both for longitude of ventilation) for ranges of values of filamentation number, F. See section 4.3 for full details of calculation.</figcaption></figure>
<figure><img src="image_26.png" /></figure>

at all during its transport into the subsurface, the squared, normalized tracer gradient would be approxi-
mately $e^{-20}$, which is close to the maximum of the distribution for $F < 15$. As the mapping is produced at
greater values of $F$, there becomes an increasing presence of large gradients. The PDFs converge toward a
maximum at $e^0$, a value that corresponds to adjacent particles originating from across the width of the
basin. As for the power spectra, the PDFs are bounded by the resolution of the mapping, since $\Delta x$ has a
lower bound of around 1 m, meaning that $\nabla C$ cannot exceed the width of the basin.

Both power spectra of tracer variance and PDFs of tracer gradients evolve in the expected manner based
on our understanding of the filamentation number, *F*. Future work may seek to elucidate a deeper relation-
ship between *F* and the properties of the tracer mapping. Such analysis could involve adopting an even
higher spatial resolution for the mapping, other metrics of the tracer distribution such as its fractal nature
*[Pierrehumbert, 1994; Neufeld et al., 2000]*, or including a form of diffusion along trajectories, as in *Abraham
and Bowen* [2002].

**5. Discussion and Conclusions**

The turbulent nature of ocean circulation is thought to play an important role in the ventilation process, but its impact remains difficult to quantify in both observations and numerical simulations. In this study, we propose that a nondimensional “filamentation number” characterizes the chaotic nature of ocean ventilation pathways. The filamentation number is the ratio of a ventilation timescale (the time since a fluid parcel was last in the mixed layer) and a Lagrangian strain timescale (the inverse of the mean strain rate following a water parcel from the mixed layer into the ocean subsurface). It quantifies the extent to which the filament width of a fluid parcel has been reduced since leaving the mixed layer, which, in turn, characterizes the extent to which pathways into the ocean interior are laminar or chaotic.

We have shown in an eddy-permitting ocean model that, due to the vastly longer timescale over which the
ocean is ventilated compared to that over which it is strained, the filamentation number is large everywhere
across three density surfaces in the subtropical North Atlantic thermocline. Consistent with this, mapping
surface properties onto these density surfaces reveals highly filamented structure. Consequently, adjacent
fluid parcels are found to have followed completely different pathways from the surface mixed layer, and
often to have been ventilated decades and great distances apart.

There is a significant increase in the filamentation number with depth (Figure 6), implying increasingly cha-
otic ventilation of the deeper density surfaces, a result confirmed by mapping the ventilation year of indi-
vidual trajectories (Figure 9). This is partly due to the longer ventilation times at depth, but also to the
decrease in $\tau_{\text{strain}}^l$ with depth (Figures 5d–5f), reflecting the different pathways along which each surface is
ventilated. These results are consistent with previous studies that have noted the increased importance of
mesoscale variability in the dynamics of the lower thermocline due to the fact that the strength of the time-
mean flow decreases more rapidly with depth than does the mesoscale eddy activity [Rhines and Schopp,
1991; Henning and Vallis, 2004]. In addition, the outcrop locations of the two deeper density surfaces
($\sigma_0=26.5 \text{ kgm}^{-3}$ and $\sigma_0=27 \text{ kgm}^{-3}$), within and to the north of the Azores current, means that direct sub-
surface pathways into the subtropical ocean are limited [Robbins et al., 2000]. Consequently, much of the
ventilation of these density surfaces must be achieved by eddy stirring across the current, leading to fila-
mentation of ventilated fluid parcels and chaotic trajectories.

In the ocean, mesoscale stirring will rapidly reduce the filament width of a coherent fluid patch down to
scales at which smaller-scale processes begin to dominate. At this stage, the ventilated fluid patch loses
coherence [Rhines and Schopp, 1991]. By characterizing the reduction in filament width since ventilation,
the filamentation number allows us to estimate the distance over which a fluid patch will retain coherence
during the ventilation process. Observations of tracers in the thermocline suggest that the filament struc-
ture generated by mesoscale eddies begins to be eroded when the filaments are reduced to widths of
around 2 km [Ledwell et al., 1998; Smith and Ferrari, 2009]. Thus, for a ventilated fluid patch with initial hori-
zontal scale, $\Delta x(0)$, on the order of 100 km (the upper bound of the nonlocal dynamics in our simulation), a
conservative estimate suggests that it would retain coherence only in regions where $F \le 4$. On all three den-
sity surfaces, this is within 100–300 km of the area where the surface is ventilated (Figure 6). This length
scale for the loss of coherence is shorter than previous estimates for a gyre circulation (300–1200 km)
derived from the Peclet number, which characterizes the relative impact of advection and eddy mixing on a
<figure><img src="image_27.png" /></figure>

tracer distribution [O'Dwyer et al., 2000; Rhines and Schopp, 1991]. This rapid erosion of the coherence of ventilated fluid parcels could severely limit the applicability of laminar ventilated thermocline theories, such as *Luyten et al.* [1983].

An alternative framing of the question of coherence is that only a fluid patch with a horizontal scale far exceeding 100 km will remain coherent a significant distance from the ventilation region. This is evident in the mapping of ventilation year, which displays a comparatively vast, coherent region of 2009 water on the shallowest density surface, onto which a large amount of water is ventilated each year. Deeper in the thermocline, the horizontal scale of recently ventilated water decreases and the filamentation occurs within a short distance of the outcrop. Furthermore, in mapping ventilation longitude, for which the initial horizontal scale is limited for all surfaces, filamentation is ubiquitous, as suggested by the large values of *F*. In the ocean, the same effect is seen in tracer observations. *Deser et al.* [1996] observed the persistence of a cold temperature anomaly during ventilation into the subtropical North Pacific. The anomaly occupied a region of more than 500 km, likely to far exceed the nonlocal regime in the ocean [LaCasce and Ohlmann, 2003; Koszalka et al., 2009], and thus the cold patch retained its signature during ventilation into the ocean interior. On the other hand, in a partially data-constrained numerical simulation of the Southern Ocean, *Jones et al.* [2016] noted erosion of the potential vorticity signature of subantarctic mode waters—with initial horizontal scales between 100 and 200 km—within a few years. Theories for the ventilated thermocline assume that horizontal gradients in potential vorticity are retained as water moves into the ocean interior [Luyten et al., 1983; *Williams*, 1991]. If the horizontal scale of these gradients is less than the scale over which filamentation occurs, our results imply that such an assumption would be invalid.

The 1/4° ocean model used in the present analysis does not fully resolve mesoscale eddies in the ocean, and is far from resolving submesoscale motions. At higher resolution, both the mesoscale strain and the details of the kinetic energy spectrum (the range of scales over which the dynamics are nonlocal) could be altered. This would affect both the value of *F* and the scales over which filamentation characterized by *F* occurs. Given the ongoing debate about the range over which ocean dynamics are local or nonlocal [LaCasce and Ohlmann, 2003; Lumpkin and Elipot, 2010; *McCaffrey et al.*, 2015], it is not obvious how our results would be altered in a higher resolution simulation or in the ocean. This remains a question for further investigation.

In conclusion, we have shown that a nondimensional filamentation number characterizes the chaotic nature of ocean ventilation. To our knowledge, this is the first attempt to apply a dynamical systems approach to the question of ocean ventilation. We have shown that it is a uniquely appropriate problem for such an approach, providing as it does an inherent timescale over which stirring by the turbulent flow takes place. Further work in this area could prove fruitful in providing deeper understanding of the processes establishing tracer variance (and by extension transit time distributions) in the ocean interior. Our results support the case for a central role for mesoscale strain in setting the properties of the subtropical thermocline. Furthermore, the kinematic approach clarifies the importance of the spatial scales of subduction and of the details of the ocean kinetic energy spectrum in the retention of memory of mixed layer properties during ventilation, with implications for the spreading of atmospherically-derived tracers into the ocean interior.

## Acknowledgments

GM acknowledges support from the Natural Environment Research Council, UK (Award Reference: 1366286), the Radcliffe-Graduate Scholarship from University College, Oxford, and the UK Met Office. The hindcast was carried out within the European Drakkar project, and the model outputs were kindly provided by J.M Molines and C. Talandier. We used the computational tool Ariane developed by B. Blanke and N. Grima, who also provided support and advice on its use. We acknowledge useful conversations with Joe LaCasce and Ric Williams, for which we are grateful. We are grateful to two anonymous reviewers, whose constructive feedback greatly improved the manuscript. Model output from the Drakkar hindcast used in this study can be obtained upon request to CL (camille.lique@ifremer.fr). Numerical codes used for analysis are available upon request to GM.

## References

* Abernathey, R. P., and D. Ferreira (2015), Southern Ocean isopycnal mixing and ventilation changes driven by winds, *Geophys. Res. Lett.*, **42**, 10,357–10,365, doi:10.1002/2015GL066238.

* Abraham, E. R., and M. M. Bowen (2002), Chaotic stirring by a mesoscale surface-ocean flow, *Chaos: Interdisciplinary J. Nonlinear Sci.*, **12**(2), 373–381, doi:10.1063/1.1481615.

* Antonsen, T. M., Z. Fan, E. Ott, and E. Garcia-Lopez (1996), The role of chaotic orbits in the determination of power spectra of passive scalars, *Phys. Fluids*, **8**(1996), 3094, doi:10.1063/1.869083.

* Barnier, B., et al. (2006), Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddy-permitting resolution, *Ocean Dyn.*, **56**(5–6), 543–567, doi:10.1007/s10236-006-0082-1.

* Barrier, N., J. Deshayes, A.-M. Treguier, and C. Cassou (2015), Heat budget in the North Atlantic subpolar gyre: Impacts of atmospheric weather regimes on the 1995 warming event, *Prog. Oceanogr.*, **130**, 75–90, doi:10.1016/j.pcean.2014.10.001.

* Bashmachnikov, I., F. Neves, T. Calheiros, and X. Carton (2015), Properties and pathways of Mediterranean water eddies in the Atlantic, *Prog. Oceanogr.*, **137**, 149–172, doi:10.1016/j.pcean.2015.06.001.

* Batchelor, G. K. (1959), Small-scale variation of convected quantities like temperature in turbulent fluid. Part 1. General discussion and the case of small conductivity, *J. Fluid Mech.*, **5**(1), 113–133, doi:10.1017/S0022112059000106.

* Bennett, A. F. (1984), Relative dispersion: Local and nonlocal dynamics, *J. Atmos. Sci.*, **41**(11), 1881–1886, doi:10.1175/1520-0469(1984)041<1881:RDLAND>2.0. CO;2.
<figure><img src="image_28.png" /></figure>

Bettencourt, J. H., C. Lopez, E. Hernandez-Garcia, I. Montes, J. Sudre, B. Dewitte, A. Paulmier, and V. Garcon (2015), Boundaries of the Peruvian oxygen minimum zone shaped by coherent mesoscale dynamics, *Nat. Geosci.*, **8**, 937–940, doi:10.1038/NGEO2570.

Blanke, B., and S. Raynaud (1997), Kinematics of the Pacific equatorial undercurrent: An Eulerian and Lagrangian approach from GCM results, *J. Phys. Oceanogr.*, **27**(6), 1038–1053, doi:10.1175/1520-0485(1997)027 < 1038:KOTPEU>2.0.CO;2.

Booth, J., and I. Kamenkovich (2008), Isolating the role of mesoscale eddies in mixing of a passive tracer in an eddy resolving model, *J. Geo-phys. Res.*, **113**, C05021, doi:10.1029/2007JC004510.

Brodeau, L., B. Barnier, A. M. Treguier, T. Penduff, and S. Gulev (2010), An ERA40-based atmospheric forcing for global ocean circulation models, *Ocean Modell.*, **31**(3–4), 88–104, doi:10.1016/j.oocemod.2009.10.005.

Callies, J., and R. Ferrari (2013), Interpreting energy and tracer spectra of upper-ocean turbulence in the submesoscale range (1?200 km), *J. Phys. Oceanog.*, **43**(11), 2456–2474, doi:10.1175/JPO-D-13-063.1.

Deser, C., M. A. Alexander, and M. S. Timlin (1996), Upper ocean thermal variations in the North Pacific during 1970–1991, *J. Clim.*, **9**, 1840–1855.

Döös, K. (1995), Inter–[O]cean exchange of water masses, *J. Geophys. Res.*, **100**(C7), 13,499–13,514.

Ferrari, R., and C. Wunsch (2009), Ocean circulation kinetic energy: Reservoirs, sources, and sinks, *Annu. Rev. Fluid Mech.*, **41**, 253–282, doi: 10.1146/annurev.fluid.40.111406. 102139.

Follows, M. J., R. G. Williams, and J. C. Marshall (1996), The solubility pump of carbon in the subtropical gyre of the North Atlantic, *J. Mar. Res.*, **54**(4), 605–630, doi:10.1357/0022240963213682.

Follows, M. J., T. Ito, and J. Marotzke (2002), The wind-driven, subtropical gyres and the solubility pump of CO2, *Global Biogeochem. Cycles*, **16**(4), 1113, doi:10.1029/2001GB001786.

Grégorio, S., T. Penduff, G. Sérazin, J.-M. Molines, B. Barnier, and J. Hirschi (2015), Intrinsic variability of the Atlantic Meridional overturning circulation at interannual-to-multidecadal time scales, *J. Phys. Oceanogr.*, **45**(7), 1929–1946, doi:10.1175/JPO-D-14-0163.1.

Haine, T. W. N., and T. M. Hall (2002), A generalized transport theory: Water-mass composition and age, *J. Phys. Oceanogr.*, **32**(6), 1932–1946, doi:10.1175/1520-0485(2002)032 < 1932:AGTTWM>2.0.CO;2.

Hall, T. M., T. W. N. Haine, D. W. Waugh, M. Holzer, F. Terenzi, and D. A. LeBel (2007), Ventilation rates estimated from tracers in the presence of mixing, *J. Phys. Oceanogr.*, **37**(11), 2599–2611, doi:10.1175/2006JPO3471.1.

Haller, G. (2015), Lagrangian coherent structures, *Annu. Rev. Fluid Mech.*, **47**(1), 137–162, doi:10.1146/annurev-fluid-010313-141322.

Henning, C. C., and G. K. Vallis (2004), The effects of mesoscale eddies on the main subtropical thermocline, *J. Phys. Oceanogr.*, **34**, 2428–2443, doi:10.1175/JPO2639.1.

Holzer, M., F. W. Primeau, W. M. Smethie, and S. Khatiwala (2010), Where and how long ago was water in the western North Atlantic ventilated? Maximum entropy inversions of bottle data from WOCE line A20, *J. Geophys. Res.*, **115**, C07005, doi:10.1029/2009JC005750.

Jenkins, W. J. (1988), The use of anthropogenic tritium and helium-3 to study subtropical gyre ventilation and circulation, *Philos. Trans. Ser. A Math. Phys. Eng. Sci.*, **325**(1583), 43–61.

Jones, D. C., A. J. S. Meijers, E. Shuckburgh, J.-B. Sallée, P. Haynes, E. K. McAufield, and M. R. Mazloff (2016), How does Subantarctic Mode Water ventilate the Southern Hemisphere subtropics?, *J. Geophys. Res. Oceans*, **121**, 6558–6582, doi:10.1002/2016JC011680.

Kamenkovich, I., Z. D. Garraffo, R. Pennel, and R. A. Fine (2017), Importance of mesoscale eddies and mean circulation in ventilation of the Southern Ocean, *J. Geophys. Res. Oceans*, **122**, 2724–2741, doi:10.1002/2016JC012292.

Karstensen, J., L. Stramma, and M. Visbeck (2008), Oxygen minimum zones in the eastern tropical Atlantic and Pacific oceans, *Prog. Oceanogr.*, **77**(4), 331–350, doi:10.1016/j.pocean. 2007.05.009.

Keating, S. R., K. S. Smith, and P. R. Kramer (2011), Diagnosing lateral mixing in the upper ocean with virtual tracers: Spatial and temporal resolution dependence, *J. Phys. Oceanogr.*, **41**, 1512–1534, doi:10.1175/2011JPO4580.1.

Koch-Larrouy, A., R. Morrow, T. Penduff, and M. Juza (2010), Origin and mechanism of Subantarctic Mode Water formation and transformation in the Southern Indian Ocean, *Ocean Dyn.*, **60**(3), 563–583, doi:10.1007/s10236-010-0276-4.

Koszalka, I. M., J. H. LaCasce, and K. A. Orvik (2009), Relative dispersion in the Nordic Seas, *J. Mar. Res.*, **67**, 411–433, doi:10.1357/002224009790741102.

LaCasce, J. H. (2002), On turbulence and normal modes in a basin, *J. Mar. Res.*, **60**(3), 431–460, doi:10.1357/002224002762231160.

LaCasce, J. H., and C. Ohlmann (2003), Relative dispersion at the surface of the Gulf of Mexico, *J. Mar. Res.*, **61**(3), 285–312, doi:10.1357/002224003322201205.

Ledwell, J. R., A. J. Watson, and C. S. Law (1998), Mixing of a tracer in the pycnocline, *J. Geophys. Res.*, **103**(C10), 21,499–21,529, doi:10.1029/98JC01738.

Lumpkin, R., and S. Elipot (2010), Surface drifter pair spreading in the North Atlantic, *J. Geophys. Res.*, **115**, C12017, doi:10.1029/2010JC006338.

Luyten, J. R., J. Pedlosky, and H. M. Stommel (1983), The ventilated thermocline, *J. Phys. Oceanogr.*, **13**(2), 292–309, doi:10.1175/1520-0485(1983)013 < 0292:TVT>2.0.CO;2.

Madec, G. (2008), NEMO Ocean Engine, Note du Pôle de modélisation, No 27, ISSN No 1288-1619, Institut Pierre-Simon Laplace (IPSL), France.

McCaffrey, K., B. Fox-Kemper, and G. Forget (2015), Estimates of ocean macroturbulence: Structure function and spectral slope from Argo profiling floats, *J. Phys. Oceanogr.*, **45**(7), 1773–1793, doi:10.1175/JPO-D-14-0023.1.

Mecking, S., M. J. Warner, C. E. Greene, S. L. Hautala, and R. E. Sonnerup (2004), Influence of mixing on CFC uptake and CFC ages in the North Pacific thermocline, *J. Geophys. Res.*, **109**, C07014, doi:10.1029/2003JC001988.

Mukiibi, D., G. Badin, and N. Serra (2016), Three dimensional chaotic advection by mixed layer baroclinic instabilities, *J. Phys. Oceanogr.*, **46**, 1509–1529, doi:10.1175/JPO-D-15-0121.1.

Neufeld, Z., P. H. Haynes, and G. Picard (2000), The effect of forcing on the spatial structure and spectra of chaotically advected passive scalars, *Phys. Fluids*, **12**(10), 2506–2513, doi:10.1063/1.1289504.

O'Dwyer, J., R. G. Williams, J. H. LaCasce, and K. G. Speer (2000), Does the potential vorticity distribution constrain the spreading of floats in the North Atlantic?, *J. Phys. Oceanogr.*, **30**, 721–732, doi:10.1175/1520-0485(2000)030 < 0721:DTPVDC>2.0.CO;2.

Pierrehumbert, R. T. (1994), Tracer microstructure in the large-eddy dominated regime, *Chaos Solitons Fractals*, **4**(6), 1091–1110, doi: 10.1016/0960-0779(94)90139-2.

Primeau, F. W. (2005), Characterizing transport between the surface mixed layer and the ocean interior with a forward and adjoint global ocean transport model, *J. Phys. Oceanogr.*, **35**, 545–564, doi:10.1175/JPO2699.1.

Primeau, F. W., and M. Holzer (2006), The ocean's memory of the atmosphere: Residence-time and ventilation-rate distributions of water masses, *J. Phys. Oceanogr.*, **36**(7), 1439–1456, doi:10.1175/JPO2919.1.

Rhines, P. B., and R. Schopp (1991), The wind-driven circulation: Quasi-geostrophic simulations and theory for nonsymmetric winds, *J. Phys. Oceanogr.*, **21**(9), 1438–1469, doi:10.1175/1520-0485(1991)021 < 1438:TWDCQG>2.0.CO;2.
<figure><img src="image_29.png" /></figure>

Rhines, P. B., and W. R. Young (1982), A theory of wind-driven circulation. *I. Mid-ocean gyres*, *J. Mar. Res.*, **40**, 559–596.

Robbins, P. E., J. Price, W. Owens, and W. J. Jenkins (2000), The importance of lateral diffusion for the ventilation of the lower thermocline in the sub-tropical North Atlantic, *J. Phys. Oceanogr.*, **30**, 67–89, doi:10.1175/1520-0485(2000)030 < 0067:TIOLDF>2.0.CO;2.

Rypina, I. I., I. Kamenkovich, P. Berloff, and L. J. Pratt (2012), Eddy-induced particle dispersion in the near-surface North Atlantic, *J. Phys. Oceanogr.*, **42**(12), 2206–2228, doi:10.1175/JPO-D-11-0191.1.

Sabine, C. L., et al. (2004), The Oceanic Sink for Anthropogenic CO2, *Science*, **305**(2004), 367–371, doi:10.1126/science.1097403.

Samelson, R. M., and G. K. Vallis (1997), Large-scale circulation with small diapycnal diffusion: The two-thermocline limit, *J. Mar. Res.*, **55**, 223–275, doi:10.1357/0022240973224382.

Shao, A. E., S. Mecking, L. Thompson, and R. E. Sonnerup (2016), Evaluating the use of 1-D transit time distributions to infer the mean state and variability of oceanic ventilation, *J. Geophys. Res. Oceans*, **121**, 6650–6670, doi:10.1002/2015JC010796.

Smith, K. S., and R. Ferrari (2009), The production and dissipation of compensated thermohaline variance by mesoscale stirring, *J. Phys. Oceanogr.*, **39**(10), 2477–2501, doi:10.1175/2009JPO4103.1.

Strogatz, S. H. (1994), *Nonlinear Dynamics and Chaos*, 2nd ed., 513 pp., Westview Press, Boulder, Colo.

Thomas, M. D., A.-M. Tréguier, B. Blanke, J. Deshayes, and A. Voldoire (2015), A Lagrangian Method to Isolate the Impacts of Mixed Layer Subduction on the Meridional Overturning Circulation in a Numerical Model, *J. Clim.*, **28**, 7503–7517, doi:10.1175/JCLI-D-14-00631.1.

Tzella, A., and P. H. Haynes (2010), Smooth and filamental structures in chaotically advected chemical fields, *Phys. Rev. E*, **81**(016322), 1–15, doi:10.1103/PhysRevE.81.016322.

Waugh, D. W., and E. R. Abraham (2008), Stirring in the global surface ocean, *Geophys. Res. Lett.*, **35**, L20605, doi:10.1029/2008GL035526.

Waugh, D. W., T. W. N. Haine, and T. M. Hall (2004), Transport times and anthropogenic carbon in the subpolar North Atlantic Ocean, *Deep Sea Res., Part I*, **51**(11), 1475–1491, doi:10.1016/j.dsr.2004.06.011.

Waugh, D. W., T. M. Hall, B. I. McNeil, R. Key, and R. J. Matear (2006), Anthropogenic CO2 in the oceans estimated using transit time distributions, *Tellus, Ser. B*, **58**(5), 376–389, doi:10.1111/j.1600-0889.2006.00222.x.

Welch, P. D. (1967), The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms, *IEEE Trans. Audio Electroacoust.*, **AU-15**, 70–73.

Wiggins, S. (2005), The dynamical systems approach to Lagrangian transport in oceanic flows, *Annu. Rev. Fluid Mech.*, **37**(1), 295–328, doi: 10.1146/annurev.fluid.37.061903.175815.

Williams, R. G. (1991), The role of the mixed layer in setting the potential vorticity of the main thermocline, *J. Phys. Oceanogr.*, **21**, 1803–1814.

Williams, R. G., J. C. Marshall, and M. A. Spall (1995), Does Stommel's Mixed Layer "Demon" Work?, *J. Phys. Oceanogr.*, **25**(12), 3089–3102, doi:10.1175/1520-0485(1995)025 < 3089:DSMLW>2.0.CO;2.

Williams, R. G., V. Roussenov, D. Smith, and M. S. Lozier (2014), Decadal evolution of ocean thermal anomalies in the north Atlantic: The effects of ekman, overturning, and horizontal transport, *J. Clim.*, **27**(2), 698–719, doi:10.1175/JCLI-D-12-00234.1.