2. Motivation: Lagrangian parcels advected by eddy-permitting state estimates

In this section we discuss the FPTD resulting from Lagrangian analysis of virtual parcels advected with eddy-permitting time varying 3-D velocity fields from two state estimates: the Southern Ocean state estimate (SOSE, Swierczek et al. Reference Swierczek, Mazloff, Morzfeld and Russell2021) and the Mediterranean Sea reanalysis (MED REA) (Simoncelli et al. Reference Simoncelli, Fratianni, Pinardi, Grandi, Drudi, Oddo and Srdjan2014). Both state estimates assimilate in-situ vertical profiles of temperature and salinity from ship-based instruments and autonomous Argo floats, in addition to sea-surface height and sea-surface temperature from satellites. The data are assimilated into ocean–sea-ice models with closed budgets of mass, momentum, temperature, salinity and sea ice, by adjusting initial conditions, and also boundary conditions for SOSE, to minimise the observations–model misfit. Both eddy-permitting models advect tracers and momentum with no explicit horizontal diffusion or viscosity: there is explicit biharmonic horizontal diffusion and viscosity with coefficients $-1\times 10^{8}$ m $^4$ s⁻¹ for SOSE and $-3\times 10^{9}$ m $^4$ s⁻¹ to remove the accumulation of gradients at the grid scale. But the unavoidable and unknown implicit diffusion due to the finite difference advection schemes probably dominates on scales larger than the grid size. The background vertical diffusion is kept to a minimum ( $1.2\times 10^{-6}$ m $^2$ s⁻¹ in MED REA and $1\times 10^{-5}$ m $^2$ s⁻¹ in SOSE). Various dynamically controlled enhancements of vertical mixing are included to avoid unstable stratification.

The details of the SOSE Lagrangian analysis are given in Rousselet, Cessi & Mazloff (Reference Rousselet, Cessi and Mazloff2023). Briefly, about eighty thousand parcels are seeded in the Atlantic sector near the equator in the upper branch of the AMOC ( $6.7$ $^\circ$ S above 1200 m), and then followed backwards in time to either the Drake Passage (66 $^\circ$ W at any depth) or to a section between the tip of South Africa and Antarctica (22 $^\circ$ E at any depth). The SOSE velocity field is resolved at $1/6$ $^\circ$ in the horizontal and with 52 vertical levels. Because 79 % of parcels originate from the 22 $^\circ$ E section, we focus on the FPTD for this group, which has the best statistics.

The details of the MED REA Lagrangian analysis are given in Vecchioni et al. (Reference Vecchioni, Cessi, Pinardi, Rousselet and Trotta2023). Briefly, about four million parcels are seeded at the Strait of Gibraltar ( $5.25$ $^\circ$ W at any depth) and then followed backwards in time in the Western Mediterranean Sea to either the Gulf of Lions, or the Tyrrhenian Sea (both at 42.31 $^\circ$ N at any depth), or the Strait of Sicily (36.94 $^\circ$ N at any depth). The MED REA velocity field is resolved at $1/16$ $^\circ$ in the horizontal and with 72 vertical levels. Because 86 % of parcels originate from the Gulf of Lions section, we focus on the FPTD for this group, which resolves even the tail of the distribution.

In both the SOSE and MED REA, each parcel is assigned a transport at release, by dividing the total transport across the initial grid face by the number of parcels released at that face. Because the flow is incompressible, each parcel’s initial transport is conserved following its trajectory.

Figure 1 shows the ensemble-averaged vertically integrated streamfunctions as a function of latitude and longitude for the two types of trajectories considered: from the 22 $^\circ$ E section to the equatorial region in the South Atlantic as represented in the SOSE (panel b); from the Strait of Gibraltar to the Gulf of Lions as represented in the MED REA (panel a). The ensemble-averaged streamfunction is constructed by summing the components of the horizontal transport carried by each parcel over all depths and occurrences in latitude–longitude bins and then summing the resulting east–west transport in latitude or the north–south transport in longitude (Blanke et al. Reference Blanke, Arhan, Madec and Roche1999). This procedure produces a two-dimensional streamfunction in the latitude–longitude plane. Thus, some 3-D spiralling motions are shown as closed contours in the Mediterranean Sea, even though each parcel path is open in three dimensions, connecting the entry and exit sections.

The two streamfunctions are qualitatively different in that the typical paths of the AMOC in the South Atlantic are direct and squeezed into a western boundary current along the South American coast (figure 1 b). Conversely, the Mediterranean paths are tortuous, punctuated by spiralling recirculations around several permanent gyres (figure 1 a). An example of the portion of a single trajectory shown in figure 2 illustrates the 3-D spiralling motion on the sub-basin-scale gyres, larger than the 50 km diameter typical of mesoscale eddies in the Western Mediterranean.

Figure 1. Ensemble-averaged and vertically integrated horizontal streamfunctions of Lagrangian parcels connecting two sections (displayed forward in time). (a) Paths from the Gulf of Lions (42.31 $^\circ$ N section in red) to the Strait of Gibraltar (5.25 $^\circ$ W section in red) as represented in MED REA in units of Sv (black contours, but contours greater than 1.35 Sv are omitted). (b) Paths from a section at 22 $^\circ$ E (connecting the tip of South Africa to Antarctica, thick line in black) to the equatorial region in the South Atlantic (6.7 $^\circ$ S section denoted by the green thick line) as represented in SOSE in units of Sverdrups ( $10^6\, \rm m^3\,s^{- 1}$ , Sv). Because there are only a few trajectories that take the retroflecting paths south of 50 $^\circ$ S, the streamfunction shows the individual trajectories.

Figure 2. Portion of a trajectory in the south Western Mediterranean, 280 days long, sampled every eight days. The arrows indicate the direction between two contiguous samples, the colour of the filled symbols indicates depth according to the colourbar on the right-hand side. The initial position is marked by a filled diamond and the final position by a filled star. Intermediate positions are marked by filled circles.

The corresponding Lagrangian (i.e. transport-weighted) FPTDs are displayed in figures 3 and 4 for the SOSE (bottom) and MED REA (top). Both distributions are compared with inverse Gaussians by fitting the first two moments (Hall et al. Reference Hall, Haine and Waugh2002; Waugh et al. Reference Waugh, Hall and Haine2003; Peacock & Maltrud Reference Peacock and Maltrud2006; Chouksey et al. Reference Chouksey, Griesel, Eden and Steinfeldt2022). Specifically, the first and second moment are calculated with the transport-weighted maximum likelihood estimator (Shao et al. Reference Shao, Mecking, Thompson and Sonnerup2016), i.e.

(2.1) \begin{equation} t_a=\frac {\sum _n\! q_n t_n}{\sum _n\! q_n}, \qquad t_d^{-1}= \frac {\sum _n\! q_n\big(t_n^{-1}-t_a^{-1}\big)}{2\sum _n\! q_n}, \end{equation}

where $t_n$ and $q_n$ are the arrival time and transport of each parcel, respectively. For the MED REA parcels, the inverse Gaussian works well with $t_{a}\approx 8$ years and $t_{d}\approx 13$ years. In this case, the effective dispersion time is comparable to the advection time across the domain.

In contrast, the fit to an inverse Gaussian (red curve in figures 3 b and 4 b) for SOSE is unsatisfactory. A better fit is provided by the Levy distribution (1.3) (blue curve in figures 3 b and 4 b). The important differences between the two theoretical distributions, most apparent in figure 4 b, is the lack of an exponential tail in (1.3) and the shift of the distribution by $t_{a}$ . We ascribe the difference in the MED REA versus SOSE FPTDs not to different statistics of mesoscale eddies in the two regions, but to the trapping of trajectories by permanent gyres in the MED REA, which are absent in the SOSE velocity field for the domain considered. In § 5 we show that the Levy distribution represents the FPTD of parcels advected by a narrow current with diffusion across the current into a region much wider than the current’s width, as appropriate for the configuration in SOSE. Conversely, in § 4 the effect of trapping by the 3-D spiralling motion associated with permanent gyres is illustrated considering stochastic processes advected by a 3-D velocity of intermediate complexity plus noise. The details of the 3-D velocity field are given in § 3.

Figure 3. Lagrangian histograms of first-passage transit times weighted by the parcels transport for the Gulf of Lions to Strait of Gibraltar trajectories in MED-RAN (panel a, black line) and for the 22 $^\circ$ E section south of the tip of Africa to the equator in the Atlantic sector in SOSE (panel b, black line). The Lagrangian histograms are compared with inverse Gaussians of the form (1.1), with the two parameters $t_a$ and $t_d$ fitted using the maximum likelihood estimator (MLE) as in (2.1) (red lines) or with a Levy distribution as in (1.3) (blue line). The parameters of the distributions, defined in the text, are reported in the titles of each panel: red numbers are for the inverse Gaussian and blue numbers are for the Levy distribution. The two parameters of the Levy distribution are obtained by the least-square fit.

Figure 4. Same as figure 3, except in logarithmic scales.

3. A 3-D steady velocity field with ocean gyres and a boundary current throughflow

To understand the interplay of permanent gyres and throughflow boundary currents on Lagrangian (or tracers) transit-time distributions, we construct a 3-D steady velocity field that includes both processes. The velocity field is the solution of the linearised, steady, primitive equations for a barotropic fluid, driven by wind stress and damped by bottom friction. We then add to this steady 3-D velocity field stochastic fluctuations to model the effects of random eddies. The resulting velocity field contains the multiscale processes that advect Lagrangian parcels in more complex contexts, such as the SOSE and MED-RAN. Each process can be studied in isolation or in different combinations to understand the resulting FPTDs.

3.1. Velocities in a barotropic layer

To illustrate the kinematics of Lagrangian trajectories advected by oceanic large-scale flows, we consider the simplest model of wind-driven gyres.

The steady, linear dynamics of a single layer of constant density, $\rho$ , and constant depth, $H$ , is governed by

(3.1) \begin{align} -\!f v &= -\frac {p_x}{\rho } -\textit{Ru} +\frac {\tau _z(y,z)}{\rho }, \end{align}

(3.2) \begin{align} \textit{fu} &= -\frac {p_y}{\rho } -\textit{R}v, \end{align}

(3.3) \begin{align} p_z &= -\rho g, \end{align}

(3.4) \begin{align} u_x&+v_y+w_z=0, \end{align}

where $f=f_0+\beta y$ , $R$ and $\rho$ are constant. The boundary conditions are that the normal velocity vanishes on the top, bottom, eastern and western boundaries, while the northern and southern boundaries are open to allow for a throughflow, $V$ , i.e.

(3.5) \begin{align} u&=0\quad {\mbox{at}}\, x=0, L_x; \end{align}

(3.6) \begin{align} v&=V(x) \quad{\mbox{at}}\, y=0, L_y; \end{align}

(3.7) \begin{align} w&=0 \quad {\mbox{at}}\, z=0, -H. \end{align}

The wind-stress forcing, $\tau (y,z)$ , is applied as a body force in the $x$ direction, defined in (3.12), and decays exponentially away from $z=0$ , representing the stress due to the wind forcing localised near the surface, while omitting the details of the Ekman spiral.

In the case of wind stress in the zonal direction only and independent of $x$ , the approximate solution for the velocity is easily derived in the limit $R\ll \beta L_x$ , and we can write the velocity field as (Rousselet & Cessi Reference Rousselet and Cessi2022)

(3.8) \begin{equation} u = -\psi _y, \qquad v = \psi _x -\phi _z, \qquad w=\phi _y, \end{equation}

where $\psi$ is the horizontal streamfunction, given by

(3.9) \begin{equation} \psi =- \frac {{\tau }_{y}(y,0)L_X}{\rho \beta H}\left ( \frac {x}{L_x}-1+ {\rm e}^{-x/\epsilon }\right) -\psi _0{\rm e}^{- x/\epsilon }, \end{equation}

where

(3.10) \begin{equation} \epsilon \equiv \frac {R}{\beta } \end{equation}

is the width of the western boundary current. The term proportional to $\psi _0$ is an exact homogeneous solution of the governing equations, as long as the throughflow velocity, $V(x)$ at $y=0,L_y$ , has the appropriate exponential form, i.e. $V(x)$ is proportional to the last term on the right-hand side of (3.9). This boundary current throughflow is a representation of the upper branch of the mid-depth overturning circulation. The strength of the gyres is controlled by the term proportional to $\tau _y$ in (3.9).

The Ekman transport is associated with an overturning streamfunction, $\phi$ , given by

(3.11) \begin{equation} \phi = \frac {\tau (y,z)}{\rho f}- \frac {\tau (y,0)}{\rho f }\left (1+\frac {z}{H}\right)\! . \end{equation}

Because of the Ekman flow near the surface and its geostrophically balanced return below the Ekman layer, $v$ has a horizontally divergent component, $\phi _z$ , which is associated with a vertical velocity. While the horizontal streamfunction is proportional to the wind-stress curl, the meridional overturning streamfunction is proportional to the wind stress. In this way, horizontal gyres are coupled at every depth by the Ekman transport and its return. Because $\psi$ and $\phi$ are approximately in quadrature in the latitudinal direction, regardless of where a parcel starts its journey, it will eventually sample all gyres and move in 3-D spiralling motions, while drifting through the domain because of the boundary current forced at the meridional boundaries, $V$ , proportional to $\psi _0$ .

The effect of permanent gyres is illustrated by choosing a wind stress of the form

(3.12) \begin{equation} \tau (y,z)\equiv \tau _0 \frac {y}{L_y}\left (\frac {y}{L_y}-1\right)\sin \left (2 \pi \frac {y}{L_y}\right){\rm e}^{z/\delta }, \end{equation}

which vanishes and has zero derivative at $y=0,L_y$ . The throughflow velocity at the boundaries $y=0,L_y$ is given by $V(x)=\epsilon ^{-1}\psi _0\exp (-x/\epsilon)$ , where $\epsilon$ is defined in (3.10).

4. Lagrangian first passage times advected by 3-D flow and stochastic noise

To illustrate the trapping role of gyres, we follow Lagrangian trajectories advected by the 3-D incompressible velocity in (3.8), (3.9) and (3.11), with an additive stochastic component proportional to $({\rm d}\xi, {\rm d}\eta, {\rm d}\zeta)$ . In other words we solve the discrete equations

(4.1) \begin{equation} {\rm d}x =u\, {\rm d}t +{\rm d}\xi, \qquad {\rm d} y = v\, {\rm d}t +{\rm d}\eta, \qquad {\rm d}z=w\, {\rm d}t + {\rm d}\zeta . \end{equation}

The stochastic component $({\rm d}\xi, {\rm d}\eta, {\rm d}\zeta)$ is picked randomly from a normal distribution with zero mean and standard deviation $(\sqrt {2{\rm d}t D}, \sqrt {2{\rm d}t D}, \sqrt {2{\rm d}t k_v})$ . Additionally, by enforcing reflection of parcels on the east and west boundaries, every trajectory bounces off the solid walls. However, diffusive flux, in addition to the prescribed advection, is allowed on the entry and exit boundaries at $y=0,L_y$ .

Table 1. Reference values of the parameters used in the barotropic model and in the Lagrangian trajectory calculations. The additional parameter $\psi _0$ is varied throughout the Lagrangian calculations.

Figure 5 shows $\psi$ and $\phi$ resulting from the wind stress in (3.12). The parameter values used are listed in table 1, and the net north–south transport, $\psi _0 H$ , is 18 Sv (1 Sverdrup – Sv = $10^6\, \rm m^3\,s^{- 1}$ ). The maximum northward western boundary current transport associated with the gyres only (i.e. without the throughflow transport proportional to $\psi _0$ ) is also 18 Sv, while the maximum southward western boundary current transport is $-12$ Sv (gyres only). In this example, the total transport in the western boundary current is everywhere northward. However, for $\psi _0$ less than a critical value that controls the strength of the gyres, the total meridional velocity on the western boundary would change sign, i.e. it would be southward in the tropical and subpolar latitudes. This threshold, whose parameter dependence is derived in § 4.2, has implications for the advective paths across the meridional length of the domain.

Figure 5. Ekman overturning streamfunction, $L_x\phi$ , as a function of latitude and depth (panel a – contour interval $3\times 10^{-3}$ Sv, negative values are dashed). Vertically integrated horizontal streamfunction, $H\psi$ , as a function of longitude and latitude (panel b – contour interval $2$ Sv, zero value at $x=L_x$ ). The expressions for $\phi$ and $\psi$ are given in (3.11) and (3.9), respectively. The parameter values are given in table 1.

As in the Lagrangian analysis using state estimates in § 2, the first passage time is defined as the first time a parcel starting at $y=0$ crosses $y=L_y$ , at any longitude, $x$ , or depth, $z$ . Parcels are initialised at $y(t=0)=0$ uniformly in $z$ . The initial release for the $n$ th parcel along $x$ is according to $x_n(t=0) =\beta /R\log (1/(1-(n-1)/N))$ ( $N$ is the number of parcels per unit depth), to ensure constant meridional transport among releases along the $x$ direction. This flux-weighted strategy ensures the equivalence of Lagrangian parcels and tracer release (Yang & Cessi Reference Yang and Cessi2024, Appendix A).

Because the steady velocity associated with the gyral circulation is three dimensional, in general the deterministic system associated with (4.1) is chaotic (Dombre et al. Reference Dombre, Frisch, Greene, Hénon, Mehr and Soward1986). Without stochastic noise, the trajectories are sensitive to initial conditions and generically fill the whole volume of the domain, although we cannot exclude small islands of periodicity (Rousselet & Cessi Reference Rousselet and Cessi2022). Thus, trajectories entering the gyres will generically sample the whole domain, distributing the throughflow velocity associated with $\psi _0$ over a wide swath of longitudes and depths, even without the stochastic noise component.

In the following we focus on the parameter values given in table 1. The Lagrangian trajectories can be examined in various limits.

4.1. Advection by boundary current plus noise, no gyres

In the absence of gyral flow, i.e. $\tau _0=0$ , the flow reduces to

(4.2) \begin{align} \psi &= -\psi _0{\rm e}^{-x/\epsilon }, \end{align}

(4.3) \begin{align} \phi &=0 . \end{align}

In this case, we have Brownian motion with an $x$ -dependent drift in the $y$ direction confined near the boundary at $x=0$ , constant in $y$ and $z$ . It is tempting to neglect the diffusion by stochastic noise and estimate the FPTD considering the statistics of parcels with advection only. The resulting FPTD is derived in Appendix A and is given by

(4.4) \begin{equation} \mbox{FPTD}(t) = {\mathcal {H}}(t-t_{a})\, \frac {t_{a}}{t^2}, \end{equation}

where $t_a\equiv \epsilon L_y/\psi _0$ . Note the $t^{-2}$ tail compared with the $t^{-3/2}$ tail for the Levy1S distribution (1.3): the latter has a `heavier’ tail because diffusion outside the western boundary layer traps the parcels for long times, and thus delays their meridional progress.

Figure 6. First passage time distributions (normalised to unity area integral) for the Lagrangian trajectories of (4.1), with the velocity field given in (3.9), for the parameter values in table 1, except $\tau _0=0$ , e.g. there are no gyres (light blue histograms with 100 bins). Results are shown for (a) $\psi _0=1.29\times 10^4$ m $^2$ s⁻¹ and (b) $\psi _0=2.57\times 10^4$ m $^2$ s⁻¹. Also shown are the theoretical predictions (5.11) (blue lines) and (4.4) (black dashed lines), using $V\equiv \psi _0/\epsilon$ , i.e. the maximum value of the boundary current velocity. The value of the parameters of the Levy1S distribution, $t_s$ and $t_a$ in years, are given in the figure.

Figure 6 shows the FPTDs for direct numerical solutions of (4.1) for two values of $\psi _0$ (light blue histograms), with all the other parameter values given in table 1. The FPTD (4.4), shown as a dashed black line, is not a good approximation for the tail of the distribution. In contrast, the Levy distribution (1.3) (blue line) is an excellent fit, with $t_a$ and $t_s$ given by

(4.5) \begin{equation} t_a=L_y\epsilon /\psi _0, \qquad t_s= D L_y^2/(\psi _0)^2, \end{equation}

using the values tabulated in table 1.The model problem developed in § 5, with a constant pipe flow diffusively leaking into a quiescent region is an appropriate idealisation for this configuration, which allows us to determine the dependence of $t_a$ and $t_s$ on the parameter values. In both examples shown in figure 6, $t_a\gt t_s$ , diffusive effects cannot be neglected and the purely advective regime is not a good fit.

4.2. Gyres plus throughflow, no noise

Another interesting limit is when there are gyres plus a throughflow boundary current, but there is no stochastic noise, i.e. $\tau _0\neq 0$ , and $D=k_v=0$ . The corresponding FPTD is shown in figure 7. The distribution is multimodal, with a slowly decaying tail. There is a large peak near $t_{a}=L_y\epsilon /\psi _0=0.9$ years in the case of the stronger boundary current (figure 7 b), which is absent for the weaker boundary current (figure 7 a). The disappearance of the peak at short transit times is controlled by the parameter

(4.6) \begin{equation} \lambda \equiv \frac {\min [-\tau _y(y,z=0)] L_x}{\rho \beta H \psi _0}\, . \end{equation}

When $|\lambda |\lt 1$ , the velocity at the western boundary is positive at every latitude. In the absence of noise, this guarantees some direct paths from $y=0$ to $y=L_y$ not held up by the gyres. This direct path is responsible for the peak at transit times $\approx 1$ year in figure 7 b. Conversely, when $|\lambda |\gt 1$ , every path from $y=0$ to $y=L_y$ is held up in the gyres, substantially lengthening the minimum transit time, as seen in figure 7 a. For the parameter values in table 1, the critical value is $\psi _0=1.85\times 10^{4}$ m $^2$ s⁻¹. In both regimes the tail of the distributions does not become exponential and decays with power laws (not shown). The qualitative difference in the trajectories depending on $\lambda$ is illustrated by examples in the two regimes in figure 8. For $|\lambda |\lt 1$ (figure 8 a), parcels starting very close to the western boundary at $x=0$ flow straight through the domain, without excursions into the gyres, with short transit times (filled circles in figure 8 a). Parcels starting further east spiral around gyres before exiting, leading to long transit times (filled diamonds in figure 8 a). For $|\lambda |\gt 1$ (figure 8 b), every trajectory experiences tortuous 3-D excursions into the gyres, leading to long transit times even for parcels starting very close to $x=0$ (e.g. filled circles in figure 8 b).

Multimodal distributions in first passage times for Lagrangian parcels have been documented in numerical ocean models of the global circulation and interpreted as arrivals of different water masses or from well-separated locations (Peacock & Maltrud Reference Peacock and Maltrud2006; Chouksey et al. Reference Chouksey, Griesel, Eden and Steinfeldt2022). Here the multimodality occurs because of the different paths taken by parcels as they avoid or get trapped into the gyres, while all starting nearby, all concentrated in the western boundary current. Once in the gyres, escapes occur in spiralling motions whose vertical component is associated with the overturning streamfunction, $\phi$ , and thus, quite slow.

Figure 7. First passage time distributions (normalised to unity area integral) for the Lagrangian trajectories of (4.1), with the velocity field given in (3.9) and (3.11), for the parameter values in table 1, except $D=k_v=0$ , e.g. there is no stochastic noise. Results are shown for (a) $H\psi _0=11$ Sv, which corresponds to $|\lambda |=1.18$ ; (b) $H\psi _0=18$ Sv, which corresponds to $|\lambda |=0.72$ . Here T10, T50 and T90 indicate the 10-percentile, 50-percentile (median) and 90-percentile transit times of the distributions (in years).

Figure 8. Pairs of Lagrangian trajectories all starting at $y=0$ , $z=-350\,\rm m$ for different $x(t=0)$ . The initial positions are marked by a green dot at $x=5$ km and a magenta diamond at $x=200$ km, and all the other parameters are as in table 1, except that $D=\kappa _v=0$ . Results are shown for (a) $H\psi _0= 18$ Sv, which corresponds to $|\lambda |=0.72$ ; (b) $H\psi _0= 11$ Sv, which corresponds to $|\lambda |=1.18$ . The exit times at $y=7000$ km are given in the textbox. The depth of the parcel is colour-coded according to the colourbar to the right of the panels.

Figure 9. First passage time distributions (normalised to unity area integral) for the Lagrangian trajectories of (4.1), with the velocity field given in (3.9) and (3.11). All parameter values are given in table 1. Results are shown for (a) $H\psi _0=11$ Sv, which corresponds to $|\lambda |=1.18$ ; (b) $H\psi _0=18$ Sv, which corresponds to $|\lambda |=0.72$ . The red curves show inverse Gaussians in (1.1) with parameters fitted by (2.1), given by $t_a$ and $t_d$ . Here T10, T50 and T90 indicate the 10-percentile, 50-percentile (median) and 90-percentile transit times of the distributions (in years).

4.3. Gyres, boundary currents and noise

When stochastic noise is added to the gyral and throughflow velocities, using the reference values in table 1, the qualitative nature of the FPTD depends on the parameter values. For $|\lambda |\lt 1$ , the deterministic velocity at the western boundary is positive at every latitude: a peak at transit times of the order of the fastest advective time is thus present even with noise (figure 9 b). When $|\lambda |\gt 1$ , all advective paths have to go through the gyral circulation. In this regime, the one-dimensional inverse Gaussian distribution is a good fit (figure 9 a, red curve). Gyres enhance the dispersion in the meridional direction beyond the diffusivity values typically expected for mesoscale eddies. In the example in figure 9 a the stochastic noise variance corresponds to a diffusivity $D=$ 500 m $^2$ s⁻¹ (cf. table 1), but the fitted $t_d=386$ years corresponds to an effective diffusivity of 4000 m $^2$ s⁻¹. Note that the dispersive enhancement by gyres is mediated by the explicit diffusivity $D$ (cf. figures 6 a, 7 a and 9 a). In the regime $|\lambda |\gt 1$ , gyres homogenise the effective throughflow meridional velocity so that the inverse Gaussian distribution, which assumes a constant throughflow, is a good approximation to the FPTD.

Trapping in the gyres delays the transit times relatively to the advective time associated with the throughflow and, with explicit diffusion, gyres collectively enhance the dispersion. For example, in figure 9(a) the advective and diffusive time scales of the inverse Gaussian are comparable, even though the Peclet number of the flow, as measured by $\max (\psi)/\!D$ , is much greater than one (in this case $\max (\psi)/\!D=25$ ).

In the range of parameters $|\lambda |\gt 1$ , where there are no direct exit paths, all parcels go through the gyres and the FPTD is well represented by the inverse Gaussian in (1.1). We expect the advection time, $t_{a}$ , to be determined by the zonally averaged throughflow, i.e. the gyres act to homogenise across the basin the net meridional transport associated with the boundary current. In this case, the theoretical estimate is $t_{a}=L_y/\bar {v}$ , where $\bar {v}$ is the zonally averaged meridional velocity. Figure 10 compares this prediction (solid line) with the estimates best fitting the FPTD from the parcel simulations (filled circles). The prediction using the homogenised velocity is qualitative, but there are quantitative discrepancies, presumably because the homogenisation is not complete and not uniform in $y$ . In these simulations all parameters are held fixed, while $\psi _0$ is varied. Also shown are the best fitted $t_{d}$ , which vary little as $\psi _0$ varies, indicating that the effective macroscopic diffusion depends on the stochastic noise and the gyral strength only.

Figure 10. The advective time, $t_{a}$ ( $\bullet$ symbol), and diffusive time, $t_{d}$ ( ${\textbf{x}}$ symbol), that best fit the invGauss FPTD in (1.1) as a function of the zonally average meridional velocity $\bar v$ . Only values of $\bar {v}$ with $|\lambda |\gt 1$ are considered. The solid curve is the theoretical prediction based on $L_y/\bar v$ . The diffusive time associated with the stochastic noise alone would be $L_y^2/\!D=3100$ years.

5. The Levy distribution for flow in a `leaky’ pipe

As a model of the behaviour observed in the SOSE (figures 1 b and 3 b), and in the intermediate complexity models of § 4.1, we construct the FPTD for parcels advected by a narrow meridional boundary current and diffused across the current. Inspired by the `hairy-pipe’ model of Young (Reference Young1988), we consider a two-dimensional problem in the $x-y$ plane such that north–south advection ( $y$ direction) is confined to a narrow region coupled with continuous diffusive `leakage’ out of the boundary current in the east–west ( $x$ ) direction. In the following, we show that this configuration results in FPTDs of the Levy form (1.3).

The probability distribution, $P$ , of parcels at a point $(x,y)$ and time $t$ is governed by

(5.1) \begin{equation} P_t+v P_y=D(P_{xx}+P_{yy}), \end{equation}

where $D$ is the diffusivity and

(5.2) \begin{equation} v(x)=\begin{cases} V &{\mbox{for}} \, -\!\epsilon \le x\le 0,\\ 0 &{\mbox{for}} \, 0\lt x \le L_x . \end{cases} \end{equation}

The velocity is a top hat confined to the narrow region $-\epsilon \le x\le 0$ , which represents the width of the boundary current, and constant within this region. Considering the domain to be $-\epsilon \le x\le L_x$ and $-\infty \le y \le L_y$ , the initial condition is

(5.3) \begin{equation} P(x,y,t=0)=\begin{cases} \epsilon ^{-1}\delta (y) &{\mbox{for}} \, -\!\epsilon \le x\le 0,\\ 0 &{\mbox{for}} \, 0\lt x \le L_x . \end{cases} \end{equation}

This flux-weighted release ensures the equivalence of Lagrangian parcels and tracer concentration distribution and is applied in the Lagrangian computations summarised in § 4 (Yang & Cessi Reference Yang and Cessi2024, Appendix A). Additional boundary conditions are $P_x=0$ at $x=-\epsilon, L_x$ and $P\to 0$ at $y=\pm \infty$ . The goal is to find the FPTDs for parcels (or tracer) starting at $y=0$ and arriving at $y=L_y$ .

Within the narrow boundary layer, $-\epsilon \le x \le 0$ , we can consider the concentration of parcels to be well mixed in the $x$ direction, while being advected by the constant meridional velocity $V\gt 0$ . Thus, for $-\epsilon \le x \le 0$ , $P=B(y,t)$ . Because the boundary layer is narrow, $\epsilon \ll L_x, L_y$ , we can also neglect diffusion in the $y$ -direction. We can determine the evolution equation for $B$ by integrating (5.1) in $x$ over $-\epsilon \le x \le 0$ to find that

(5.4) \begin{equation} \epsilon (B_t + \textit{VB}_y)= \textit{DP}_x(x=\epsilon, y,t) . \end{equation}

The assumption of a well-mixed concentration in the boundary layer holds as long as the advection is slow compared with the diffusion, i.e. $V/L_y \ll D/\epsilon ^2$ .

For $x\ge 0$ , there is no advection and just diffusion of tracer. Anticipating a solution that decays rapidly away from $x=0$ , we neglect diffusion in the $y$ direction and $P$ satisfies

(5.5) \begin{equation} P_t = \textit{DP}_{xx}. \end{equation}

Despite the zero initial condition in the region $x\ge 0$ , there is a concentration $P$ introduced by the diffusive flux leaking out of the boundary layer, i.e.

(5.6) \begin{equation} P(x=0,y,t)=B(y,t). \end{equation}

Integrating (5.5) in $x$ for $0\le x\le L_x$ we find that

(5.7) \begin{equation} \int _0^{L_x} P_{t}\, {\rm d}x=-\textit{DP}_{x}(x=0, y,t) . \end{equation}

Summing (5.4) and (5.7) we obtain the evolution of the integral of $P$ over the entire $x$ domain, governed by

(5.8) \begin{equation} \epsilon ( B_{ t}+\textit{VB}_{y}) + \int _0^{L_x} P_{t}\, {\rm d}x=0. \end{equation}

Because we neglect diffusion in the $y$ direction, and given the initial condition localised at $y=0$ , no parcels are expected to be in the $y\lt 0$ part of the domain (assuming that $V\gt 0$ ). Therefore, $P$ and $B$ vanish for $y\le 0$ . The distribution of first passage time (FPTD) at the generic point $y$ starting from the initial release at $y=0$ is given by (minus) the change in the survival probability (Cox & Miller Reference Cox and Miller1965, Chapters 5.4, 5.7), i.e.

(5.9) \begin{equation} \mbox{FPTD}(y,t)=- \int _0^y\!\! {\rm d}y^{\prime} \left [ \epsilon\! B_{ t}(y^{\prime},t) + \int _0^{L_x}\!\!\!\! P_{t} (x,y^{\prime},t) \, {\rm d}x \right] \! . \end{equation}

The right-hand side of (5.9) can be obtained by integrating (5.8) in $y$ , and using the boundary condition $B=0$ at $y=0$ , to obtain

(5.10) \begin{equation} \mbox{FPTD}(y,t)= \epsilon\! \textit{VB}(y,t) . \end{equation}

In order to find $B(y,t)$ , we need to solve (5.4) in the region $x\lt 0$ with $P_{x}(x=0, y,t)$ determined by solving the interior diffusion (5.5) for $ x\gt 0$ . The solution for the FPTD, denoted with `Levy1S’, can be found using the Laplace transform in the limit $L_x/\epsilon \to \infty$ (Appendix B) and is given by

(5.11) \begin{equation} \mbox{Levy1S}(y=L_y,t;t_a;t_s) ={\mathcal {H}}(t-t_{a})\sqrt {\frac {t_s}{4\pi (t-t_{a})^3}} \exp \left [-\frac {t_s}{4(t-t_{a})}\right], \end{equation}

where ${\mathcal {H}}(x)$ is the Heaviside function. This is the same expression as (1.3), with explicit expressions for the two time scales as a function of the physical parameters of the model. The advective time, $t_a$ , and the scale parameter, $t_s$ , are defined as

(5.12) \begin{equation} t_{a}\equiv L_y/V, \qquad t_s\equiv \frac {D L_y^2}{(\epsilon V)^2}=t_{a}^2/t_{D}, \end{equation}

where $t_{D}\equiv \epsilon ^2/\!D$ is the diffusive time across the boundary current. Because diffusion in the $y$ direction is neglected, no parcel can exit the domain before the advective time $t_{a}$ . After this time, the distribution changes on the time scale $t_s$ . This expression for the FPTD is a (one-sided) Levy distribution, with the shift parameter given by $t_{a}$ and the scale parameter given by $t_s$ . Because of the slow $t^{-3/2}$ decay of Levy1S for large times, all moments of the distribution diverge, so we cannot calculate the average first passage time. For times $t\gt t_{a}$ , (5.11) coincides with the FPTD of Brownian motion without drift, as appropriate in the limit $V/L_y \ll D/\epsilon ^2$ .

Figure 6 shows that the distribution Levy1S (blue line) describes well the FPTD of parcels advected by a combination of stochastic noise plus a western boundary current that decays exponentially from the boundary, using $V=\psi _0/\epsilon$ , i.e. the maximum velocity of the current. The fit (5.11) works well, even though $V/L_y=(1.8, 3.7)\times 10^{-8}\,{\mbox{s}}^{-1}$ is not much smaller than $D/\epsilon ^2= 5 \times 10^{-8}\,{\mbox{s}}^{-1}$ , indicating that this simplified model has a range of validity beyond the limit in which it was derived.

Figure 6 compares the numerically simulated distribution (light blue histogram) to the purely advective solution (4.4) (black dashed line) and the Levy’s distribution (5.11) (blue line). Even though the advective solution (4.4) uses the exact expression for the profile of the boundary current, it is a poorer agreement than the Levy’s distribution (5.11) using a constant velocity for the boundary current. This is because the parcels slowly leaking diffusively in the quiescent interior are important for the weakly decaying tail of the FTPD.

6. A two-peaked distribution with two diffusively coupled pipes

Motivated by the results in § 4.3, we further expand the two-dimensional model developed in § 5, by adding a weak meridional velocity outside the western boundary current region. This `outer’ meridional velocity represents the slow net meridional movement of parcels caught in the recirculating gyres. The rapid velocity in the western boundary current represents the fast meridional movement of parcels that can proceed without recirculating around the gyres. In the three-gyre model of § 4.3 these fast paths only exist for $|\lambda |\lt 1$ .

With these assumptions, the probability distribution, $P$ , of parcels at a point $(x,y)$ and time $t$ is governed by

(6.1) \begin{equation} P_t+v P_y=D^x\!P_{xx}+D^y\!P_{yy}, \end{equation}

where

(6.2) \begin{equation} v(x)=\begin{cases} V_0 &{\mbox{for}} \, -\!\epsilon \le x\le 0,\\ V_1 &{\mbox{for}} \, 0\lt x \le L_x . \end{cases} \end{equation}

The domain is $-\epsilon \le x\le L_x$ and $-\infty \le y\le L_y$ , where $\epsilon$ is the width of the fast throughflow western boundary current. We consider the diffusivity to be anistropic, because the tortuous paths induced by the gyres enhance the diffusivity in the $y$ direction. To account for parcels that are deviated by the gyral circulation very close to the entry point at $y=0$ (cf. figure 8 b), we modify the initial conditions to

(6.3) \begin{equation} P= \begin{cases} \delta (y)&{\mbox{for}} \, -\!\epsilon \le x\le 0,\\ \mu \delta (y)&{\mbox{for}} \, 0\lt x \le L_x . \end{cases} \end{equation}

The parameter $\mu$ is the strength of the release applied in the large interior region. Considering that the gyral circulation deviates an $O(1)$ fraction of the parcels injected at $y=0$ into the interior, we consider $\mu =O(\epsilon /L_X)$ , so that the total transport in the boundary current and in the outer region are comparable. Additional boundary conditions are $P_x=0$ at $x=-\epsilon, L_x$ and $P\to 0$ as $y\to -\infty$ .

Within the narrow boundary layer, $-\epsilon \le x \le 0$ , the concentration of parcels is well mixed in the $x$ direction, while being advected by the constant velocity meridional velocity $V_0\gt 0$ . Because the boundary layer is narrow, $\epsilon \ll L_x, L_y$ , we can also neglect diffusion in the $y$ -direction. Thus, for $-\epsilon \le x \le 0$ , $P=B(y,t)$ . We can determine the evolution equation for $B$ by integrating (5.1) in $x$ over the boundary layer region to find that

(6.4) \begin{equation} \epsilon (B_t + V_0B_y)= D^x\! P_x(x=0, y,t) . \end{equation}

For $x\ge 0$ , there is weak advection and diffusion of tracer and $P$ satisfies

(6.5) \begin{equation} P_t + V_1\! P_y = D^x\!P_{xx}+D^y\!P_{yy}. \end{equation}

The concentration must be continuous at $x=0$ , so we have

(6.6) \begin{equation} P(x=0,y,t)=B(y,t) \, . \end{equation}

The details of the approximate solution can be found in Appendix C. The final expression for the FPTD at the point $L_y$ is

(6.7) \begin{align} &\text{FPTD}(L_y,t)={\mathcal {H}}(t-t_0){\mathcal {H}}(t_1-t)\frac {\sqrt {t_0}[r(t_1-t) +2r\sqrt {t_0/t_D}(t-t_0)]}{2\sqrt {\pi (1-r)(t-t_{0})^3}}\exp \big (- \xi ^2 \big)\nonumber\\ &\quad +\mu {\mathcal {H}}(t-t_0){\mathcal {H}}(t_1-t)\left [\frac {\sqrt {t_0}(t-t_0)}{\sqrt {\pi (1-r)(t-t_{0})^3}}\exp (- \xi ^2) +\frac { r\sqrt {t_D}}{\sqrt {t_0}(1-r)}\mbox{erfc}(\xi)\right] \nonumber\\ &\quad +\mu \sqrt {\frac {t_y t_0}{4\pi t^3}}\exp \left [- \frac {t_y(t-t_1)^2}{4 t_1^2t}\right]\left [\sqrt { t_x}\mbox{erf}\left (\frac {t_x}{2\sqrt {t}}\right) -2\sqrt {\frac {t}{\pi }} \right]\! .\end{align}

We have defined

(6.8) \begin{align} \xi \equiv \frac { r(t-t_1)}{2\sqrt {t_D(t-t_0)}},\qquad r\equiv V_1/V_0, \qquad t_0 &\equiv L_y/V_0, \qquad t_1\equiv L_y/V_1, \end{align}

(6.9) \begin{align} t_{D}\equiv \epsilon ^2/D^x, \qquad t_y &\equiv L_y^2/ D^y, \qquad t_x\equiv L_x^2/D^x . \end{align}

Here $r$ is the ratio of the meridional velocities, $t_0$ is the advective time along the fast boundary current and $t_1$ is the advective time in the slow interior region. These two advective time scales control the two peaks of the distribution. In (6.7), $t_{D}$ is the rapid diffusion across the boundary layer, $t_y$ and $t_x$ are the slow diffusive times along the full latitudinal and longitudinal extents of the domain, respectively. The diffusive time scales $t_{D}$ and $t_y$ control the widths of the two peaks. The first term on the right-hand side of (6.7) is a modification to the expression for the LevyS1 (5.11) due the slow advection outside the boundary current, proportional to $r\ll 1$ . The long time $t_1$ would allow the tail of the distribution to decay exponentially for large times. However, without meridional diffusion (neglected to obtain this term), all parcels exit the domain after a time $t_1$ and none reach the absorbing boundary before the fast advection time $t_0$ . Thus, the exponential decay for $t\gg t_1$ is not sampled by this term. The second term on the right-hand side of (6.7) accounts for the backscatter towards $x\gt 0$ of parcels originally initialised in the interior then diffused back into $x\lt 0$ . The last term on the right-hand side of (6.7) is the inverse Gaussian arising from the advection by $V_1$ in the interior and the diffusion in the latitudinal direction, described by the parameter $t_y$ . Diffusion in the longitudinal direction, controlled by the parameter $t_x$ , plays a minor role: the contribution to the structure in the longitudinal direction is essentially constant in the parameter range of interest.

An example of the FPTD for a reasonable choice of parameters is shown in figure 11. We do not attempt a detailed fit to the FPTD obtained by direct numerical integrations of the stochastic system (4.1). The important point is that there are two peaks in the distribution: the peak at short times, similar to the Levy1S distribution, is associated with the advection along streamlines in the western boundary current that do not recirculate in the interior gyres, e.g. the trajectory with filled circles in figure 8 a. The peak at later times is associated with the streamlines that recirculate around gyres and are advected by a slower meridional velocity distributed across the whole domain, while experiencing enhanced diffusion in the $y$ -direction associated with the stirring effect of the large-scale gyres. The analysis in § 4.3 indicates that the gyres increase the effective meridional diffusivity by a factor larger than one (but less than ten). In figure 11, $D^y$ is increased by a factor of about five relative to $D^x$ to illustrate this process.

Figure 11. First passage time distribution (6.7). The values of the dimensional parameters are $\epsilon =100$ km, $\mu =0.0429$ , $D^x=500$ m $^2$ s⁻¹, $D^y=2390$ m $^2$ s⁻¹, $V_0=0.129$ m s⁻¹, $V_1=0.0055$ m s⁻¹, $L_x=L_y=7\times 10^{3}$ km. These values correspond to $t_0=1.73$ years, $t_1=40.28$ years, $t_x=650$ years, $t_y=3100$ years. The values of the non-dimensional parameters introduced in Appendix C are reported in the legend.

The distribution exhibits two peaks: one at $t=t_0$ , associated with the first term on the right-hand side of (6.7), and one at $t=t_1$ , associated with the inverse Gaussian on the third line of (6.7). The discontinuity at $t=t_1$ is associated with the corresponding Heaviside function in the FTPD in (6.7), and would be smoothed out by the addition of diffusion in the $y$ direction in the boundary layer region. When all the streamlines recirculate in the gyres, the width of the boundary region vanishes and we can set $t_D= 0$ . Then, the only term that survives is the term on the third line of (6.7) proportional to the inverse Gaussian distribution.