Ocean transit times: from basin to planetary scales

TL;DR

The paper explores how ocean currents and eddies affect the time it takes for water to travel across large ocean regions.

Contribution

The study reveals how multiscale processes, including throughflow currents and eddy dispersion, shape Lagrangian transit times in the ocean.

Findings

Mesoscale eddies increase dispersion and widen transit-time distributions.

Ocean gyres cause tortuous paths and multimodal transit-time distributions.

Trapping within gyres leads to complex transit-time patterns.

Abstract

Lagrangian transit times on basin to planetary scales are controlled by the interplay of multiscale processes. The primary advective time scale is set by throughflow currents, such as interhemispheric western boundary currents. Dispersion by mesoscale eddies introduces fluctuations that erase memory and enhance dispersion, widening the transit-time distribution. The tortuous paths of Lagrangian parcels, particularly within ocean gyres, significantly enhance dispersion beyond the levels attributed to mesoscale eddies alone. Additionally, trapping by ocean gyres leads to multimodal distributions of Lagrangian transit times. These processes are illustrated in three complementary contexts: eddy-permitting ocean state estimates, simplified spatially extended three-dimensional flows and diffusively coupled two-dimensional pipe models.

Figures

Full text

Introduction

Information about the time scales and pathways between regions in the ocean can be obtained by either tracking the advection-diffusion of tracers or following Lagrangian trajectories. Both methods are used in observational campaigns, using transient tracers (Fine 2011) or Lagrangian floats (Bower et al. 2019). The observational studies have been complemented by numerical studies, where tracers (Haine & Hall 2002; Holzer & Primeau 2007; DeVries & Primeau 2011; Shah et al. 2017) or parcels are advected and diffused using synthetic velocity fields produced by ocean models (Döös 1995; Blanke & Raynaud 1997; Speich, Blanke & Cai 2007; Rousselet, Cessi & Forget 2021).

The statistics of tracers and Lagrangian parcels, represented by either tracer concentrations or probability distributions of the parcel’s density in space and time, are equivalent as long as the boundary and initial conditions are consistent (Shah et al. 2017). The Lagrangian method allows efficient computations in the limit where the diffusion time is much longer than the advection time, a regime that is difficult to achieve with Eulerian methods. In the following, we focus on the first passage time statistics derived from Lagrangian analysis. First passage time is defined as the time when a stochastic walker or a diffusing parcel first reaches a given site starting from another given site (Redner 2001).

In the oceanic context, the dispersion by mesoscale eddies is a well-studied feature in Lagrangian analysis, which can be quantified in terms of the diffusivity, $[eqn]$ , (Davis 1991a,* b *). The effect of eddies might be considered to be analogous to the diffusion induced by stochastic noise in probability space in the Fokker–Planck (or Kolmogorov) equation describing the probability of parcels being at a point in space at a given time, when advected by a random velocity. One implicit assumption is that eddies’ have smaller scales, in space or time, than the mean’ component of the velocity. Additional assumptions are that eddies transport tracers down the mean gradient, at least in some direction, e.g. along isopycnals.

Oceanic flows on the planetary scale have multiscale motions including throughflow currents that are typically confined to a boundary, but span the full meridional extent of basins across two hemispheres, e.g. the western boundary currents associated with the Atlantic meridional overturning circulation (AMOC). On the basin scales are large-scale gyres, which also include narrow boundary currents on the same scale of the throughflow and of mesoscale eddies. Gyres, in particular, trap trajectories into spiralling three-dimensional (3-D) motions (Rousselet et al. 2021; Rousselet & Cessi 2022; Vecchioni et al. 2023). Mesoscale eddies diffuse tracers efficiently along isopycnals, delaying the progression of parcels and tracers inside boundary currents. The net effect of gyres and eddies is to broaden the transit-time distributions of tracers and parcels travelling meridionally. In this work we focus on the upper ocean circulation, i.e. the flow in the top kilometre of the water column, which is the best-observed region of the ocean. We attempt to tease apart the contribution of throughflows, gyres and eddies to the transit times between sites separated in latitude.

Several works successfully fit the first passage time distribution (FPTD) of Eulerian tracers or Lagrangian parcels in the ocean as an inverse Gaussian distribution of the form (Hall, Haine & Waugh 2002; Waugh, Hall & Haine 2003; Peacock & Maltrud 2006; Chouksey et al. 2022)

$$[eqn]$$

The inverse Gaussian (1.1) arises in the probability distribution of first passage times for parcels, in one dimension, starting at $[eqn]$ and arriving at $[eqn]$ subject to a uniform drift, $[eqn]$ , and Brownian motion with zero mean and variance given by $[eqn]$ (Cox & Miller 1965, Chapter 5.7). The form (1.1) clarifies that the parameters of the inverse Gaussian can be interpreted in terms of advective, $[eqn]$ , and diffusive, $[eqn]$ , times associated with the underlying Fokker–Planck equation. The two parameters are given by

$$[eqn]$$

In (1.1) the relevant diffusive time is along the direction of the constant velocity, $[eqn]$ , while in higher dimensions, diffusion across the mean velocity is another relevant parameter.

The inverse Gaussian distribution is characterised by a small but finite probability at short times. This is because rare but large fluctuations of the stochastic component can take parcels from the starting point to the arrival target in one random step. Equivalently, for finite but non-zero diffusivity, an exponentially small but finite tracer concentration is diffused immediately along a line, before the advective propagation. At large times, the inverse Gaussian distribution decays exponentially, because the advection can efficiently flush out parcels (or tracers), despite the diffusive backscatter.

An alternative to the inverse Gaussian distribution is the (one-sided) Levy distribution, given by

$$[eqn]$$

where $[eqn]$ is the Heaviside function. As for the inverse Gaussian, the Levy distribution contains two parameters: the shift parameter (sometimes called location’ parameter), $[eqn]$ , represents a threshold time that needs to be overcome before the target site can ever be reached; the so called scale’ parameter, $[eqn]$ , determines the position of the distribution maximum relative to $[eqn]$ (Nolan 2020). Unlike the inverse Gaussian, Levy1S vanishes for short times, $[eqn]$ . Contrary to the inverse Gaussian, the Levy distribution decays slowly for large times, as $[eqn]$ . This heavy tail causes all the moments larger than zero to diverge. We show in § 5 that the distribution (1.3) is appropriate for first passage times of parcels in two or three dimensions, advected by a narrow throughflow current with minimum advection time $[eqn]$ , and subject to noise/diffusion across the throughflow current, into a reservoir much wider than the width of the throughflow current. Because of the diffusive hold-up into the reservoir, parcels’ progression is slowed down and the tail of the FPTD is heavy.

In § 2 it is shown that the two distributions, invGauss and Levy1S, arise in two separate regions of the ocean, the South Atlantic Ocean and the Western Mediterranean Sea, where different processes control the dynamics.

To our knowledge, the Levy1S distribution has not been discussed in the oceanographic context. The trapping by excursions in regions of temporary retention, such as gyres, has been studied in the context of transport of tracers (Young 1988) and porous media (Hidalgo, Neuweiler & Dentz 2021). Here, we revisit the problem in the oceanographic context, adding together throughflow boundary currents on the planetary scale, basin-scale gyres and stochastic noise (parametrising mesoscale eddies) in the advection of Lagrangian parcels, i.e. the main elements of the circulation in the top kilometre of the ocean. In § 4 it is shown that the distribution invGauss is appropriate for first passage times when, in the presence of eddies, the western boundary currents associated with gyres overcome the throughflow boundary current. In this case, the effective diffusivity in the direction of the throughflow boundary current is enhanced above the value estimated from eddy dispersion alone. Conversely, the Levy1S distribution dominates the fast time component of the transit-time distribution when the throughflow boundary currents velocity overcomes the opposing western boundary current associated with the gyral transport (also in the presence of eddies).

This study also shows that large-scale gyres lead to multimodal distributions, in addition to holding up parcels for long periods. This implies that multimodal distributions are not necessarily due to different origins of parcels or water masses, but rather to different routes of parcels with contiguous origins.

The retarding and dispersive effects of gyres and eddies is illustrated using a hierarchy of models of decreasing complexity, with the common process of Lagrangian parcels advection. The most complex models are eddy-permitting state estimates used to advect virtual parcels in basin-scale regions that accommodate three types of flows: (i) currents that flow through the open boundaries of the domain; (ii) permanent gyres on scales comparable to the domain, which recirculate parcels in 3-D spiralling motions; (iii) transient mesoscale eddies on scales much smaller than the region of interest. All three types of flows contribute to the trajectories connecting the open boundaries of the region of interest. The FPTDs in these complex 3-D, time-dependent flows and in complicated geometrical domains are discussed in § 2: they motivate the formulation and analysis of simplified models.

The intermediate complexity model uses the 3-D steady velocity solution of the linear barotropic momentum equations in a semi-enclosed domain forced by a simplified surface wind stress with three gyres plus a north–south throughflow boundary current flowing across the open boundaries. The throughflow current is a crude representation of the interhemispheric throughflow boundary current associated with the upper branch of the AMOC. The formulation of this velocity field is detailed in § 3. The total velocity field advects Lagrangian parcels, which are also subject to stochastic Brownian motion to account for time-dependent flows on smaller scales. This framework allows us to explore different regimes of relative importance of gyres versus thoughflow currents, and the role of eddies (noise) in enabling the gyral dispersive effect. The results of the stochastic Lagrangian advection is discussed in § 4.

The simplest models consider constant flows through pipes that are diffusely coupled in two dimensions. This simplification allows us to derive the explicit form of the FPTD as a function of the parameters of the model. The pipe models are discussed in §§ 5 and 6.

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. 2021) and the Mediterranean Sea reanalysis (MED REA) (Simoncelli et al. 2014). 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 $[eqn]$ m $[eqn]$ s^−1^ for SOSE and $[eqn]$ m $[eqn]$ s^−1^ 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 ( $[eqn]$ m $[eqn]$ s^−1^ in MED REA and $[eqn]$ m $[eqn]$ s^−1^ 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 (2023). Briefly, about eighty thousand parcels are seeded in the Atlantic sector near the equator in the upper branch of the AMOC ( $[eqn]$ $[eqn]$ S above 1200 m), and then followed backwards in time to either the Drake Passage (66 $[eqn]$ W at any depth) or to a section between the tip of South Africa and Antarctica (22 $[eqn]$ E at any depth). The SOSE velocity field is resolved at $[eqn]$ $[eqn]$ in the horizontal and with 52 vertical levels. Because 79 % of parcels originate from the 22 $[eqn]$ 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. (2023). Briefly, about four million parcels are seeded at the Strait of Gibraltar ( $[eqn]$ $[eqn]$ 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 $[eqn]$ N at any depth), or the Strait of Sicily (36.94 $[eqn]$ N at any depth). The MED REA velocity field is resolved at $[eqn]$ $[eqn]$ 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 $[eqn]$ 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. 1999). 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 $[eqn]$ N section in red) to the Strait of Gibraltar (5.25 $[eqn]$ 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 $[eqn]$ E (connecting the tip of South Africa to Antarctica, thick line in black) to the equatorial region in the South Atlantic (6.7 $[eqn]$ S section denoted by the green thick line) as represented in SOSE in units of Sverdrups ( $[eqn]$ , Sv). Because there are only a few trajectories that take the retroflecting paths south of 50 $[eqn]$ 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. 2002; Waugh et al. 2003; Peacock & Maltrud 2006; Chouksey et al. 2022). Specifically, the first and second moment are calculated with the transport-weighted maximum likelihood estimator (Shao et al. 2016), i.e.

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the arrival time and transport of each parcel, respectively. For the MED REA parcels, the inverse Gaussian works well with $[eqn]$ years and $[eqn]$ 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 $[eqn]$ . 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 $[eqn]$ 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 $[eqn]$ and $[eqn]$ 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.

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.

Velocities in a barotropic layer

3.1.

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, $[eqn]$ , and constant depth, $[eqn]$ , is governed by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ , $[eqn]$ and $[eqn]$ 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, $[eqn]$ , i.e.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The wind-stress forcing, $[eqn]$ , is applied as a body force in the $[eqn]$ direction, defined in (3.12), and decays exponentially away from $[eqn]$ , 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 $[eqn]$ , the approximate solution for the velocity is easily derived in the limit $[eqn]$ , and we can write the velocity field as (Rousselet & Cessi 2022)

$$[eqn]$$

where $[eqn]$ is the horizontal streamfunction, given by

$$[eqn]$$

where

$$[eqn]$$

is the width of the western boundary current. The term proportional to $[eqn]$ is an exact homogeneous solution of the governing equations, as long as the throughflow velocity, $[eqn]$ at $[eqn]$ , has the appropriate exponential form, i.e. $[eqn]$ 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 $[eqn]$ in (3.9).

The Ekman transport is associated with an overturning streamfunction, $[eqn]$ , given by

$$[eqn]$$

Because of the Ekman flow near the surface and its geostrophically balanced return below the Ekman layer, $[eqn]$ has a horizontally divergent component, $[eqn]$ , 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 $[eqn]$ and $[eqn]$ 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, $[eqn]$ , proportional to $[eqn]$ .

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

$$[eqn]$$

which vanishes and has zero derivative at $[eqn]$ . The throughflow velocity at the boundaries $[eqn]$ is given by $[eqn]$ , where $[eqn]$ is defined in (3.10).

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 $[eqn]$ . In other words we solve the discrete equations

$$[eqn]$$

The stochastic component $[eqn]$ is picked randomly from a normal distribution with zero mean and standard deviation $[eqn]$ . 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 $[eqn]$ .

Table 1.Reference values of the parameters used in the barotropic model and in the Lagrangian trajectory calculations. The additional parameter $[eqn]$ is varied throughout the Lagrangian calculations.ParameterValueNotes $[eqn]$

$[eqn]$ mBasin width $[eqn]$

$[eqn]$ mBasin length $[eqn]$

$[eqn]$ mBasin depth $[eqn]$

$[eqn]$ Bottom drag $[eqn]$

$[eqn]$ mEkman layer depth $[eqn]$ Total Coriolis parameter $[eqn]$

$[eqn]$ Coriolis parameter at $[eqn]$

$[eqn]$

$[eqn]$

$[eqn]$ parameter $[eqn]$

$[eqn]$ Wind-stress amplitude $[eqn]$

$[eqn]$ Constant density $[eqn]$

$[eqn]$ Horizontal diffusion $[eqn]$

$[eqn]$ Vertical diffusion $[eqn]$

$[eqn]$ Width of western boundary layer

Figure 5 shows $[eqn]$ and $[eqn]$ resulting from the wind stress in (3.12). The parameter values used are listed in table 1, and the net north–south transport, $[eqn]$ , is 18 Sv (1 Sverdrup – Sv = $[eqn]$ ). The maximum northward western boundary current transport associated with the gyres only (i.e. without the throughflow transport proportional to $[eqn]$ ) is also 18 Sv, while the maximum southward western boundary current transport is $[eqn]$ Sv (gyres only). In this example, the total transport in the western boundary current is everywhere northward. However, for $[eqn]$ 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, $[eqn]$ , as a function of latitude and depth (panel a – contour interval $[eqn]$ Sv, negative values are dashed). Vertically integrated horizontal streamfunction, $[eqn]$ , as a function of longitude and latitude (panel b – contour interval $[eqn]$ Sv, zero value at $[eqn]$ ). The expressions for $[eqn]$ and $[eqn]$ 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 $[eqn]$ crosses $[eqn]$ , at any longitude, $[eqn]$ , or depth, $[eqn]$ . Parcels are initialised at $[eqn]$ uniformly in $[eqn]$ . The initial release for the $[eqn]$ th parcel along $[eqn]$ is according to $[eqn]$ ( $[eqn]$ is the number of parcels per unit depth), to ensure constant meridional transport among releases along the $[eqn]$ direction. This flux-weighted strategy ensures the equivalence of Lagrangian parcels and tracer release (Yang & Cessi 2024, 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. 1986). 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 2022). Thus, trajectories entering the gyres will generically sample the whole domain, distributing the throughflow velocity associated with $[eqn]$ 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.

Advection by boundary current plus noise, no gyres

4.1.

In the absence of gyral flow, i.e. $[eqn]$ , the flow reduces to

$$[eqn]$$ $$[eqn]$$

In this case, we have Brownian motion with an $[eqn]$ -dependent drift in the $[eqn]$ direction confined near the boundary at $[eqn]$ , constant in $[eqn]$ and $[eqn]$ . 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

$$[eqn]$$

where $[eqn]$ . Note the $[eqn]$ tail compared with the $[eqn]$ 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 $[eqn]$ , e.g. there are no gyres (light blue histograms with 100 bins). Results are shown for (a) $[eqn]$ m $[eqn]$ s^−1^ and (b) $[eqn]$ m $[eqn]$ s^−1^. Also shown are the theoretical predictions (5.11) (blue lines) and (4.4) (black dashed lines), using $[eqn]$ , i.e. the maximum value of the boundary current velocity. The value of the parameters of the Levy1S distribution, $[eqn]$ and $[eqn]$ in years, are given in the figure.

Figure 6 shows the FPTDs for direct numerical solutions of (4.1) for two values of $[eqn]$ (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 $[eqn]$ and $[eqn]$ given by

$$[eqn]$$

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 $[eqn]$ and $[eqn]$ on the parameter values. In both examples shown in figure 6, $[eqn]$ , diffusive effects cannot be neglected and the purely advective regime is not a good fit.

Gyres plus throughflow, no noise

4.2.

Another interesting limit is when there are gyres plus a throughflow boundary current, but there is no stochastic noise, i.e. $[eqn]$ , and $[eqn]$ . The corresponding FPTD is shown in figure 7. The distribution is multimodal, with a slowly decaying tail. There is a large peak near $[eqn]$ 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

$$[eqn]$$

When $[eqn]$ , the velocity at the western boundary is positive at every latitude. In the absence of noise, this guarantees some direct paths from $[eqn]$ to $[eqn]$ not held up by the gyres. This direct path is responsible for the peak at transit times $[eqn]$ year in figure 7 b. Conversely, when $[eqn]$ , every path from $[eqn]$ to $[eqn]$ 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 $[eqn]$ m $[eqn]$ s^−1^. 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 $[eqn]$ is illustrated by examples in the two regimes in figure 8. For $[eqn]$ (figure 8 a), parcels starting very close to the western boundary at $[eqn]$ 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 $[eqn]$ (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 $[eqn]$ (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 2006; Chouksey et al. 2022). 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, $[eqn]$ , 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 $[eqn]$ , e.g. there is no stochastic noise. Results are shown for (a) $[eqn]$ Sv, which corresponds to $[eqn]$ ; (b) $[eqn]$ Sv, which corresponds to $[eqn]$ . 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 $[eqn]$ , $[eqn]$ for different $[eqn]$ . The initial positions are marked by a green dot at $[eqn]$ km and a magenta diamond at $[eqn]$ km, and all the other parameters are as in table 1, except that $[eqn]$ . Results are shown for (a) $[eqn]$ Sv, which corresponds to $[eqn]$ ; (b) $[eqn]$ Sv, which corresponds to $[eqn]$ . The exit times at $[eqn]$ 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) $[eqn]$ Sv, which corresponds to $[eqn]$ ; (b) $[eqn]$ Sv, which corresponds to $[eqn]$ . The red curves show inverse Gaussians in (1.1) with parameters fitted by (2.1), given by $[eqn]$ and $[eqn]$ . Here T10, T50 and T90 indicate the 10-percentile, 50-percentile (median) and 90-percentile transit times of the distributions (in years).

Gyres, boundary currents and noise

4.3.

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 $[eqn]$ , 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 $[eqn]$ , 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 $[eqn]$ 500 m $[eqn]$ s^−1^ (cf. table 1), but the fitted $[eqn]$ years corresponds to an effective diffusivity of 4000 m $[eqn]$ s^−1^. Note that the dispersive enhancement by gyres is mediated by the explicit diffusivity $[eqn]$ (cf. figures 6 a, 7 a and 9 a). In the regime $[eqn]$ , 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 $[eqn]$ , is much greater than one (in this case $[eqn]$ ).

In the range of parameters $[eqn]$ , 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, $[eqn]$ , 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 $[eqn]$ , where $[eqn]$ 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 $[eqn]$ . In these simulations all parameters are held fixed, while $[eqn]$ is varied. Also shown are the best fitted $[eqn]$ , which vary little as $[eqn]$ varies, indicating that the effective macroscopic diffusion depends on the stochastic noise and the gyral strength only.

Figure 10.The advective time, $[eqn]$ ( $[eqn]$ symbol), and diffusive time, $[eqn]$ ( $[eqn]$ symbol), that best fit the invGauss FPTD in (1.1) as a function of the zonally average meridional velocity $[eqn]$ . Only values of $[eqn]$ with $[eqn]$ are considered. The solid curve is the theoretical prediction based on $[eqn]$ . The diffusive time associated with the stochastic noise alone would be $[eqn]$ years.

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 ([1988](#ref31)), we consider a two-dimensional problem in the $[eqn]$ plane such that north–south advection ( $[eqn]$ direction) is confined to a narrow region coupled with continuous diffusive leakage’ out of the boundary current in the east–west ( $[eqn]$ ) direction. In the following, we show that this configuration results in FPTDs of the Levy form (1.3).

The probability distribution, $[eqn]$ , of parcels at a point $[eqn]$ and time $[eqn]$ is governed by

$$[eqn]$$

where $[eqn]$ is the diffusivity and

$$[eqn]$$

The velocity is a top hat confined to the narrow region $[eqn]$ , which represents the width of the boundary current, and constant within this region. Considering the domain to be $[eqn]$ and $[eqn]$ , the initial condition is

$$[eqn]$$

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 2024, Appendix A). Additional boundary conditions are $[eqn]$ at $[eqn]$ and $[eqn]$ at $[eqn]$ . The goal is to find the FPTDs for parcels (or tracer) starting at $[eqn]$ and arriving at $[eqn]$ .

Within the narrow boundary layer, $[eqn]$ , we can consider the concentration of parcels to be well mixed in the $[eqn]$ direction, while being advected by the constant meridional velocity $[eqn]$ . Thus, for $[eqn]$ , $[eqn]$ . Because the boundary layer is narrow, $[eqn]$ , we can also neglect diffusion in the $[eqn]$ -direction. We can determine the evolution equation for $[eqn]$ by integrating (5.1) in $[eqn]$ over $[eqn]$ to find that

$$[eqn]$$

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. $[eqn]$ .

For $[eqn]$ , there is no advection and just diffusion of tracer. Anticipating a solution that decays rapidly away from $[eqn]$ , we neglect diffusion in the $[eqn]$ direction and $[eqn]$ satisfies

$$[eqn]$$

Despite the zero initial condition in the region $[eqn]$ , there is a concentration $[eqn]$ introduced by the diffusive flux leaking out of the boundary layer, i.e.

$$[eqn]$$

Integrating (5.5) in $[eqn]$ for $[eqn]$ we find that

$$[eqn]$$

Summing (5.4) and (5.7) we obtain the evolution of the integral of $[eqn]$ over the entire $[eqn]$ domain, governed by

$$[eqn]$$

Because we neglect diffusion in the $[eqn]$ direction, and given the initial condition localised at $[eqn]$ , no parcels are expected to be in the $[eqn]$ part of the domain (assuming that $[eqn]$ ). Therefore, $[eqn]$ and $[eqn]$ vanish for $[eqn]$ . The distribution of first passage time (FPTD) at the generic point $[eqn]$ starting from the initial release at $[eqn]$ is given by (minus) the change in the survival probability (Cox & Miller 1965, Chapters 5.4, 5.7), i.e.

$$[eqn]$$

The right-hand side of (5.9) can be obtained by integrating (5.8) in $[eqn]$ , and using the boundary condition $[eqn]$ at $[eqn]$ , to obtain

$$[eqn]$$

In order to find $[eqn]$ , we need to solve (5.4) in the region $[eqn]$ with $[eqn]$ determined by solving the interior diffusion (5.5) for $[eqn]$ . The solution for the FPTD, denoted with `Levy1S’, can be found using the Laplace transform in the limit $[eqn]$ (Appendix B) and is given by

$$[eqn]$$

where $[eqn]$ 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, $[eqn]$ , and the scale parameter, $[eqn]$ , are defined as

$$[eqn]$$

where $[eqn]$ is the diffusive time across the boundary current. Because diffusion in the $[eqn]$ direction is neglected, no parcel can exit the domain before the advective time $[eqn]$ . After this time, the distribution changes on the time scale $[eqn]$ . This expression for the FPTD is a (one-sided) Levy distribution, with the shift parameter given by $[eqn]$ and the scale parameter given by $[eqn]$ . Because of the slow $[eqn]$ decay of Levy1S for large times, all moments of the distribution diverge, so we cannot calculate the average first passage time. For times $[eqn]$ , (5.11) coincides with the FPTD of Brownian motion without drift, as appropriate in the limit $[eqn]$ .

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 $[eqn]$ , i.e. the maximum velocity of the current. The fit (5.11) works well, even though $[eqn]$ is not much smaller than $[eqn]$ , 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.

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 $[eqn]$ .

With these assumptions, the probability distribution, $[eqn]$ , of parcels at a point $[eqn]$ and time $[eqn]$ is governed by

$$[eqn]$$

where

$$[eqn]$$

The domain is $[eqn]$ and $[eqn]$ , where $[eqn]$ 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 $[eqn]$ direction. To account for parcels that are deviated by the gyral circulation very close to the entry point at $[eqn]$ (cf. figure 8 b), we modify the initial conditions to

$$[eqn]$$

The parameter $[eqn]$ is the strength of the release applied in the large interior region. Considering that the gyral circulation deviates an $[eqn]$ fraction of the parcels injected at $[eqn]$ into the interior, we consider $[eqn]$ , so that the total transport in the boundary current and in the outer region are comparable. Additional boundary conditions are $[eqn]$ at $[eqn]$ and $[eqn]$ as $[eqn]$ .

Within the narrow boundary layer, $[eqn]$ , the concentration of parcels is well mixed in the $[eqn]$ direction, while being advected by the constant velocity meridional velocity $[eqn]$ . Because the boundary layer is narrow, $[eqn]$ , we can also neglect diffusion in the $[eqn]$ -direction. Thus, for $[eqn]$ , $[eqn]$ . We can determine the evolution equation for $[eqn]$ by integrating (5.1) in $[eqn]$ over the boundary layer region to find that

$$[eqn]$$

For $[eqn]$ , there is weak advection and diffusion of tracer and $[eqn]$ satisfies

$$[eqn]$$

The concentration must be continuous at $[eqn]$ , so we have

$$[eqn]$$

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

$$[eqn]$$

We have defined

$$[eqn]$$ $$[eqn]$$

Here $[eqn]$ is the ratio of the meridional velocities, $[eqn]$ is the advective time along the fast boundary current and $[eqn]$ is the advective time in the slow interior region. These two advective time scales control the two peaks of the distribution. In (6.7), $[eqn]$ is the rapid diffusion across the boundary layer, $[eqn]$ and $[eqn]$ are the slow diffusive times along the full latitudinal and longitudinal extents of the domain, respectively. The diffusive time scales $[eqn]$ and $[eqn]$ 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 $[eqn]$ . The long time $[eqn]$ 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 $[eqn]$ and none reach the absorbing boundary before the fast advection time $[eqn]$ . Thus, the exponential decay for $[eqn]$ is not sampled by this term. The second term on the right-hand side of (6.7) accounts for the backscatter towards $[eqn]$ of parcels originally initialised in the interior then diffused back into $[eqn]$ . The last term on the right-hand side of (6.7) is the inverse Gaussian arising from the advection by $[eqn]$ in the interior and the diffusion in the latitudinal direction, described by the parameter $[eqn]$ . Diffusion in the longitudinal direction, controlled by the parameter $[eqn]$ , 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 $[eqn]$ -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, $[eqn]$ is increased by a factor of about five relative to $[eqn]$ to illustrate this process.

Figure 11.First passage time distribution (6.7). The values of the dimensional parameters are $[eqn]$ km, $[eqn]$ , $[eqn]$ m $[eqn]$ s^−1^, $[eqn]$ m $[eqn]$ s^−1^, $[eqn]$ m s^−1^, $[eqn]$ m s^−1^, $[eqn]$ km. These values correspond to $[eqn]$ years, $[eqn]$ years, $[eqn]$ years, $[eqn]$ years. The values of the non-dimensional parameters introduced in Appendix C are reported in the legend.

The distribution exhibits two peaks: one at $[eqn]$ , associated with the first term on the right-hand side of (6.7), and one at $[eqn]$ , associated with the inverse Gaussian on the third line of (6.7). The discontinuity at $[eqn]$ 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 $[eqn]$ 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 $[eqn]$ . Then, the only term that survives is the term on the third line of (6.7) proportional to the inverse Gaussian distribution.

Discussion

The different elements of the ocean general circulation, interhemispheric boundary currents, basin-scale gyres and mesoscale eddies, each contribute to a spectrum of transit times between target regions on the planetary scales. A subset of the possible rich behaviour encountered in the oceanographic context of multiscale flows is illustrated by the transit-time distributions of two disparate regions of the world ocean: the South Atlantic Ocean and the Western Mediterranean Sea.

Because of the connection with the Southern Ocean, the subtropical gyre of the southern hemisphere does not close within the Atlantic sector. Thus, upper ocean paths from the tip of South Africa to the equatorial region of the western South Atlantic are direct except for the dispersion effected by mesoscale eddies. This type of circulation is well described by a boundary current with noise. This configuration leads to FPTDs that lack exponential tails and decay with slow power laws for large times. This behaviour is in contrast with the classical paradigm of Brownian motion with constant drift, described by the inverse Gaussian, which decays exponentially for large times and has well-defined moments of the distribution.

In contrast, the Western Mediterranean Sea is connected to adjacent seas by narrow and shallow straits, which promote gyral circulations on several scales, in addition to boundary currents connecting the interior flow to the straits, especially to the Strait of Gibraltar. The net result of the interaction between these throughflow currents, the gyres and the eddies are transit-time distributions that, when integrated over the directions perpendicular to the throughflow are well described by a one-dimensional inverse Gaussian distribution.

In general, 3-D gyral circulation leads to transit-time distributions potentially quite complex, depending on the strength of the gyres relative to the throughflow currents and the mesoscale eddies. Insofar as eddies can be represented with white noise, some physical intuition can be obtained by considering Lagrangian parcels advected by steady 3-D flows plus stochastic noise. With weak or no stochastic noise, multimodal distributions are the norm, even when the origin of the parcels (or tracer) is localised to a small region. This is because the tortuosity of the paths varies among trajectories that start close together: parcels that are very close at the initial time, quickly diverge taking disparate times to cross the domain, even in the ensemble average. A fast throughflow boundary current further complicates the distribution, by adding quick paths that avoid the spiralling gyral routes.

When stochastic noise is added to throughflow boundary current plus gyres, the distribution can be unimodal or multimodal depending on the value of the parameters controlling the total velocity of the boundary current. If the boundary current changes sign along the boundary then all parcel paths recirculate within the gyres. Gyral recirculation increases the transit times by homogenising the net throughflow velocity over a broader region. At the same time, gyres increase parcel dispersion well above the typical values associated with mesoscale eddies. The net result is a FPTD that is broad and decays exponentially for large times.