The present study reconciles theoretical differences between the Lagrangian diffusivity and effective diffusivity in a transformed spatial coordinate based on the contours of a quasi-conservative tracer. In the transformed coordinate, any adiabatic stirring effect, such as shear-induced dispersion, is naturally isolated from diabatic cross-contour motions. Therefore, Lagrangian particle motions in the transformed coordinate obey a transformed zeroth-order stochastic (i.e., random walk) model with the diffusivity replaced by the effective diffusivity. Such a stochastic model becomes the theoretical foundation on which both diffusivities are exactly unified. In the absence of small-scale diffusion, particles do not disperse at all in the transformed contour coordinate. Besides, the corresponding Lagrangian autocorrelation becomes a delta function and is thus free from pronounced overshoot and negative lobe at short time lags that may be induced by either Rossby waves or mesoscale eddies; that is, particles decorrelate immediately and Lagrangian diffusivity is already asymptotic no matter how small the time lag is. The resulting instantaneous Lagrangian spreading rate is thus conceptually identical to the effective diffusivity that only measures the instantaneous irreversible mixing. In these regards, the present study provides a new look at particle dispersion in contour-based coordinates.
Estimating the lateral eddy diffusivity of a passive tracer in the ocean, revealing its temporal/spatial variabilities, and understanding its relation to flow characteristics are the key issues of subgrid-scale parameterization in numerical models. Up to now, direct and indirect estimates of lateral diffusivity may generally be categorized into three kinds.
- Eulerian diffusivity: A commonly used parameterization relates the eddy flux of a passive tracer q to its mean gradient as (e.g., Bachman et al. 2015; Fox-Kemper et al. 2013): Bachman et al. 2015). The Eulerian diffusivity (tensor) K can thus be inferred through the eddy flux and the mean gradient of q (e.g., Vallis 2006). However, extra constraints are needed to uniquely determine K components (e.g., Bachman et al. 2015; Fox-Kemper et al. 2003; Nakamura 2001), and several physically based constraints have already been proposed (e.g., Eden et al. 2009; Marshall and Shutts 1981; McDougall and McIntosh 1996; Nakamura 2001; Osborn and Cox 1972).
- Lagrangian diffusivity: Lagrangian diffusivity (tensor) is defined, after the decorrelation time scale, as the half growth rate of the displacement variance (dispersion) of Lagrangian particles (e.g., Fox-Kemper et al. 2013; LaCasce 2008): Taylor 1922) and extended to inhomogeneous oceanic flows (Davis 1987, 1991). This method requires only a large number of current-following (Lagrangian) observations and becomes popular in descriptive studies of the oceans (LaCasce 2008; Lumpkin and Pazos 2007), as well as in studies pursuing diffusivity estimates (e.g., Chiswell 2013; LaCasce et al. 2014; Peng et al. 2015; Qian et al. 2013; Sallée et al. 2008; Zhurbas and Oh 2004).
- Effective diffusivity: This definition was first given by Rhines and Young (1983), and systematically reformulated by Nakamura (1996) [also see Winters and D’Asaro (1996)]. The key concept is to adopt the area A enclosed by a contour of a conservative tracer as a spatial coordinate (contour-based coordinate). Through this coordinate transform, the advection–diffusion equation is converted into a diffusion-only equation in which the diffusivity is called effective diffusivity (Nakamura 1996; Haynes and Shuckburgh 2000): multiplied by an amplifying factor (squared ratio of equivalent length to minimum possible length for a given contour) representing the efficiency of stirring (e.g., Lee et al. 2009; Marshall et al. 2006). If tracer contours are well stirred, effective diffusivity could be greatly amplified; without stirring, however, the effective diffusivity reduces to (Nakamura 2008).
As the definitions of the above three eddy diffusivities are quite distinct, one cannot generally expect an exact correspondence from one to another. Large discrepancies among these diffusivity estimates, including both the magnitude and pattern, are also reported in the literature (e.g., Abernathey et al. 2010; Abernathey and Marshall 2013; Chiswell 2013; Klocker et al. 2012b; Lu and Speer 2010; Marshall et al. 2006; Peng et al. 2015; Riha and Eden 2011; Sallée et al. 2008; Zhurbas et al. 2014). Reconciling these discrepancies will help understand why they differ, provide a clear interpretation, and eventually improve the robustness of different estimates for parameterization. The work by Klocker et al. (2012b) was one of few studies that make great effort to reduce these discrepancies. They emphasized that only after decorrelation time scale (i.e., Lagrangian diffusivity becomes asymptotic) should the Lagrangian diffusivity be equivalent to effective diffusivity. The present study works on reconciling the Lagrangian diffusivity and the effective diffusivity, but demonstrates a new perspective of observing particle dispersion in a contour-based coordinate. This eventually leads to an exact theoretical unification of the two diffusivities, complementing the work by Klocker et al. (2012b).
The present reconciliation lies on a physical rationale that advective (reversible) and diffusive (irreversible) effects of oceanic eddies should be isolated and parameterized separately (Table 1; Gent 2011; Gent et al. 1995; Lee et al. 1997). It is known that advection acts to redistribute tracer without altering its Lagrangian properties [i.e., the material derivative (e.g., Bennett 1987)] and requires adiabatic parameterization (e.g., Gent and McWilliams 1990; Gent et al. 1995). Diffusion, on the other hand, decays tracer variance and homogenizes tracer until uniformity is reached (e.g., Marshall et al. 2006; Nakamura 2008). As a result, it can be readily parameterized by the downgradient closure (e.g., Nakamura 2001; Redi 1982).
The isolation of advection from diffusion in the context of effective diffusivity is direct, as the effective diffusivity measures only irreversible mixing (Nakamura 2008). Similar isolation also applies to the Eulerian diffusivity. It has long been known that Eulerian eddy flux can be decomposed into rotational and divergent components (e.g., Jayne and Marotzke 2002; Lau and Wallace 1979; Lee et al. 1997; Plumb 1979; Roberts and Marshall 2000), equivalent to the split of the Eulerian diffusivity tensor K into antisymmetric and symmetric parts (e.g., Bachman et al. 2015; Fox-Kemper et al. 2013; Nakamura 2001). The rotational (antisymmetric) component is equivalent to an extra advective velocity (Gent et al. 1995; Vallis 2006) whereas the divergent (symmetric) component is modeled by, in isotropic cases, a downgradient scalar diffusivity (e.g., Eden et al. 2007; Nakamura 2001).
Such a physically based isolation also has been pursued in the perspective of Lagrangian diffusivity. It is well known that stirring (more precisely, differential advection), other than diffusion, greatly amplifies Lagrangian particle dispersion. This is also known as shear dispersion (Taylor 1953; Young and Jones 1991). Therefore, particle dispersion generally contains both reversible (shear dispersion) and irreversible contributions (Table 1). Shearing flows usually result in ballistic increase of particle dispersion (e.g., Oh et al. 2000), preventing Lagrangian diffusivity from asymptotic (e.g., Bauer et al. 1998). As a result, great efforts were put on removing shear dispersion. A common way to reduce shear dispersion, in the context of single-particle statistics, is to perform a scale separation (), remove the time-mean flow , and use the eddy part to infer Lagrangian diffusivity (e.g., Bauer et al. 1998; Davis 1991; LaCasce 2008; Peng et al. 2015; Poulain 2001). However, this only removes the shear effect of . Shear effect may still reside in . Although there are studies trying to remove the shear effect within using, for example, a bicubic fitting inside a Eulerian bin (e.g., Qian et al. 2014), a spin parameter to fit the looping feature (e.g., Veneziani et al. 2004), or a principle axes decomposition of diffusivity tensor (e.g., Oh et al. 2000), shear-induced dispersion is only minimized in some sense rather than completely being removed, for example, there is a residual diffusivity (Wolfram and Ringler 2017). Relative dispersion diffusivity (two-particle statistics), although mitigating the mean flow shear (Batchelor 1952; Bennett 1987), could still be affected by the shear between the two selected particles (e.g., LaCasce and Bower 2000). This leads to a fact that Lagrangian diffusivity may also sense some parts of the reversible mixing. Since the effective diffusivity only measures irreversible mixing, this becomes a major obstacle that one has to overcome in reconciling the two diffusivities.
Based upon these physical arguments above, here we propose a novel reconciliation, in a contour-based coordinate rather than in a geographic Eulerian coordinate (Fig. 1), between Lagrangian diffusivity and effective diffusivity. In such a coordinate, particles decorrelate immediately and the Lagrangian dispersion diffusivity only measures instantaneous irreversible mixing. As a result, both two diffusivities can be exactly unified. It is worth noting that some inspiring works have already incorporated contour coordinates to investigate Lagrangian dispersion. LaCasce and Speer (1999) reported dispersion characteristics in a contour coordinate based on time-mean potential vorticity. Sallée et al. (2011) used a contour coordinate based on sea surface height to investigate the cross-stream dispersion in the Southern Ocean. In the present study, there are several nontrivial modifications from previous studies: 1) abandoning the use of periodic reentrant models, 2) using (quasi) conservative tracer contours, and 3) estimating dispersion across instantaneous contours rather than across time-mean (climatological) contours. It will be shown that these modifications provide new insights in the present theoretical reconciliation.
a. Numerical model setup
The Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al. 1997) is used to generate both tracer data and Lagrangian particle trajectories for the diffusivity calculations. The underlying flow follows the nondivergent barotropic vorticity equation on a beta plane ( m−1 s−1 and s−1 at ) in Cartesian x–y coordinates, similar to the classical setup for simulating a wind-driven double-gyre oceanic circulation. This is achieved in the primitive equation MITgcm by imposing a rigid-lid constraint. The domain is rectangular with a uniform grid of 5.5 km and closed no-slip boundary conditions, covering an area of 3080 km (x direction) × 2200 km (y direction). The bottom of the ocean is flat with a depth of 3500 m and a small linear drag. The wind stress is time independent and only varies sinusoidally in the y direction with a peak amplitude of 0.5 N m−2 at the center of the domain. A time step of 300 s is used to advance the model and the flow is spinup from rest. The Leith viscosity scheme (Leith 1996) for momentum dissipation is used as the flow dynamics falls into the regime of 2D geostrophic turbulence (Fox-Kemper and Menemenlis 2008).
The online infrastructures of MITgcm for advecting a passive tracer and Lagrangian particles are used. The passive tracer is initialized as a double-gyre like distribution (Fig. 2a) and then advected by the nondivergent flow for 1440 days (4 model years, see snapshots in Figs. 2c,d), using the second-order moment Prather scheme (Prather 1986) suggested by Hill et al. (2012). For the tracer simulation demonstrated here, an explicit diffusivity of is used. However, spurious numerical diffusivity from the Prather advection scheme is inevitable. In this case, we estimate its magnitude using the domain-integrated tracer variance equation as (e.g., Marshall et al. 2006):
where the overbar indicates domain average. Then can be obtained by observing the evolution of global-mean tracer variance and squared tracer gradient (Fig. 3a). Taken as the ratio of the above two quantities, is generally below 6 m2 s−1 (Fig. 3b), indicating that the tracer evolves quasi-adiabatically. Lagrangian particles are also deployed with a uniform interval of 4 km into the flow from the start. All these particles are tracked using the standard fourth-order Runge–Kutta method, finally resulting in more than 400 000 particles (excluding those that collide with the boundaries). Both the tracer and particle data are output daily for subsequent analyses.
The closed-domain setup here is chosen on purpose to highlight some differences from previous studies. Previous studies using effective diffusivity preferentially chose a global or hemispherical domain (e.g., Haynes and Shuckburgh 2000; Marshall et al. 2006), or zonally periodic reentrant domain (e.g., Abernathey and Marshall 2013; Ferrari and Nikurashin 2010; Klocker et al. 2012b; Shuckburgh et al. 2009a; Wolfram and Ringler 2017). Even when the zonal periodicity is not exactly satisfied, it is forced to be (e.g., Klocker et al. 2012b; Shuckburgh et al. 2009a). In this case, the zonal mean is a good approximation to the along-contour mean (e.g., Klocker et al. 2012b; Wolfram and Ringler 2017). Besides, the Lagrangian diffusivity tensor could also be simplified to a meridional component only, and this (cross-stream) component is expected to be less affected by shear dispersion induced by (e.g., LaCasce et al. 2014). The present closed-domain setup actually breaks the equivalence of the zonal mean and the along-contour mean. As will be shown later, our reconciliation here does not depend on this equivalence and is thus more general than that of Klocker et al. (2012b).
b. Coordinate transformation
The contour-based coordinate is constructed as follows (e.g., Marshall et al. 2006; Nakamura 1996, 2008; Wolfram and Ringler 2017). First, the area A enclosed by a selected contour q of the tracer field is calculated as
using the “box-counting” method (Nakamura 2008). To reduce discretization error, model-output tracer data are interpolated linearly to double the resolution before counting areas. A number of 201 equally spaced tracer contours between minimum and maximum values are used to generate the tabulated area–contour relation . Then a transformed Y coordinate position
is calculated and reinterpolated to 400 equally spaced Y grid points (same grid points as y coordinate). In practice, it is found occasionally that the contour interval is not large enough to ensure the strict monotonicity of . That is, the interval of neighboring contours is too small so that no samples fall into the contour bin. In this case, the contour bins are adjusted by looking up and interpolating the tabulated from equally spaced area bins. The procedure above eventually provides a coordinate transformation from geographic x–y coordinates to the contour-based X–Y coordinates. In the transformed coordinate, the X coordinate is trivial (no variation), and its zonal average is exactly the along-contour average in the geographic coordinates.
Tracer evolutions will also be quite different in these two coordinates, which can be appreciated from Fig. 2. As the tracer is advected quasi-adiabatically by the underlying flow, tracer contours in the x–y coordinates are stretched and elongated by the flow, or even broken into isolated small contour islands (Figs. 2c,d). However, as the area A enclosed by a tracer contour does not change with time in the adiabatic limit (Nakamura 1996), the coordinate transformation from Figs. 2c and 2d could be almost the same as Fig. 2b, indicating that the tracer evolves very slowly in the X–Y coordinate. A similar argument also applies to Lagrangian particles (the black dot in Fig. 2). In the adiabatic limit, a Lagrangian particle should remain on the same contour no matter what the shape of the contour is, and its Y position does not change with time in the contour space (also see Fig. 1).
3. Theoretical reconciliation
The theoretical arguments start from the advection–diffusion equation that the modeled passive tracer q, as demonstrated in Fig. 2, follows
where v is the nondivergent horizontal velocity and is either the spurious numerical diffusivity caused by, for example, the Eulerian advection scheme ( could vary in space but we assume this variation negligible and treat it as a constant here) or the explicit diffusivity set in the model. Dynamically consistent with Eq. (1), the well-known zeroth-order stochastic (random walk) model for the position of a single Lagrangian q molecule is
where and are the particle position and Weiner process vectors whose component is a random increment that follows the normal distribution (e.g., Berloff and McWilliams 2002; Brickman and Smith 2002; Griffa 1996; Rodean 1996). This means that Eq. (1) and Eq. (2) are exactly the same model but of different implementations, that is, Eq. (1) for the concentration of q while Eq. (2) for a single q molecule. Therefore, Eq. (1) is also called the Fokker–Planck equation for Eq. (2) (e.g., Griffa 1996).
If Eq. (3) is viewed as a one-dimensional Fokker–Plank equation and q as the positional probability for Lagrangian particles in the Y coordinate, with introspection we can write the corresponding random-walk model (e.g., Berloff and McWilliams 2002; Brickman and Smith 2002; Griffa 1996; Rodean 1996) in the Y coordinate as
Equation (4) is essentially the one-dimensional zeroth-order stochastic model that particles should obey in the transformed Y coordinate. Comparing Eq. (4) with Eq. (2), the advective term in Eq. (2), which is in the Lagrangian frame of reference, is dropped as all the advective motions have already been absorbed into the motion of the Y-coordinate surfaces (Nakamura 2008). Note that the first term on the rhs of Eq. (4), representing the effect of inhomogeneous diffusion [a.k.a. drift correction (e.g., Berloff and McWilliams 2002; Brickman and Smith 2002)], is necessary as it makes the model satisfy the well-mixed criterion (Thomson 1987). This term naturally vanishes in Eq. (2) as is assumed as a constant in space and in Eq. (4) varies on account of .
Now Eqs. (3) and (4) become the Y-coordinate counterparts respectively to Eqs. (1) and (2) in the x–y coordinates. Equation (4) also forms the theoretical basis for discussing the relation between particle dispersion and effective diffusivity. Several scenarios can be identified from this expression. First, in the absence of small-scale diffusion (), and therefore particle position remains constant in the Y coordinate. This means that in the Y-coordinate particles do not disperse at all in the absence of small-scale diffusion (Fig. 1). Second, in the case that the effective diffusivity is a constant (homogeneous scenario), Eq. (4) becomes
Squaring both sides, taking the ensemble average , and noting that , we have
where is the Lagrangian particle displacement during a time interval of . This expression directly relates particle dispersion to Nakamura’s effective diffusivity in the transformed Y coordinate and can also be written as
which states that Lagrangian single-particle diffusivity, defined as the half growth rate of dispersion (e.g., LaCasce 2008), is exactly the effective diffusivity. Therefore, dispersion in the Y coordinate is caused by small-scale diffusion () and also enhanced by adiabatic stirring, consistent with the concept of effective diffusivity (e.g., Nakamura 1996). In the third scenario, varies with space (inhomogeneous scenario, which is a common case in turbulent flows) and the mean drift (first moment) of particle displacement should be removed before obtaining . The mean drift is
given that , which indicates that the mean drift is proportional to the gradient of . Then the Lagrangian diffusivity (cloud diffusivity or relative diffusivity),
is exactly Nakamura’s effective diffusivity. Note that one could define a nondimensional Péclet number in the Y coordinate as the ratio of the mean drift and the stochastic motion:
given that follows the normal distribution of . In the limit of a large where the inhomogeneity of diffusion (gradient of stirring) dominates over the magnitude of the diffusivity itself, particles drift systematically toward regions of high diffusivity (e.g., Davis 1991), like being “advected” by a velocity of . For small cases where turbulence is relatively homogeneous, particles will spread out randomly without any systematic drift. In this case, the center of mass of the particles remains almost unchanged so that the single-particle diffusivity (or absolute diffusivity) is equivalent to the cloud diffusivity (e.g., LaCasce 2008).
It is interesting to compare the above theoretical arguments with previous findings. Klocker et al. (2012b) demonstrated in a reentrant model that the asymptotic estimates of Lagrangian diffusivity in cross-stream direction (y component) agrees well with effective diffusivity estimates. However, according to Eq. (6), it is found that only in the transformed Y coordinate are both diffusivities exactly equivalent. As their meridional dispersion contains shear dispersion induced by (i.e., their meridional displacement is caused mainly by advection of ), their estimates inevitably contain adiabatic shear contribution. This is evident in their Lagrangian autocorrelation that is characterized by a pronounced overshoot followed by a negative lobe at relatively short time lags. This feature is associated with advective processes, for example, the back and forth meandering effect of Rossby waves embedded in an eastward jet (e.g., Klocker et al. 2012a) or the looping effect of mesoscale eddies (e.g., Veneziani et al. 2004). Such meandering or looping effect should be termed as reversible mixing that are induced by (not by ), rendering the Lagrangian diffusivity at short time lags conceptually different from effective diffusivity that only measures irreversible mixing. Only after the decorrelation time scale, asymptotic Lagrangian diffusivity might be dominated by irreversible mixing, which supports the equivalence of the two. At short time lags, such equivalence cannot be expected. In addition, the reconciliation here does not depend on the reentrant domain setup because the along-contour mean in the x–y coordinates is exactly the “zonal” mean in the Y coordinate.
Sallée et al. (2011) also estimated particle dispersion in their contour coordinate based on sea surface height. As sea surface height is not exactly a conservative quantity, other adiabatic processes (e.g., divergent flow) as well as diabatic processes other than diffusion (e.g., wind forcing) will lead to its Lagrangian changes. Fortunately, they did find that particle displacement in their contour coordinate is more Gaussian than in geographical coordinates. According to Eq. (4), this is reasonable because statistical properties of particle displacement D, after removing the mean drift , should be identical to those of dW (exactly Gaussian). However, any temporal average to obtain a Gaussian probability density function (PDF) may be inadequate, as the magnitude of effective diffusivity should vary with time. Only ensemble averaging at a time instant would guarantee a Gaussian PDF since displacement amplitude is determined by instantaneous effective diffusivity. Their results confirm that sea surface height is a practical choice to minimize the advective effects of rotational eddies and meandering jets, as compared to traditional binning method.
Shuckburgh and Haynes (2003) also compared particle transport in a kinematic flow with the solution of the transformed advection–diffusion Eq. (3) [their Eq. (11)]. They found that there is a remarkably good quantitative agreement (see their Fig. 12) between particle density estimated using a simple kernel function in contour-based coordinates and that predicted by Eq. (3). They also emphasized that they had no precise theoretical understanding of why this should be the case. Our study has showed that particles in the contour-based coordinates should obey the transformed random-walk model of Eq. (4) with effective diffusivity regulating the random displacement impulse, which is precisely the theoretical foundation they were seeking. Since the corresponding Fokker–Planck equation (controlling the evolution of particle concentration) for Eq. (4) is exactly Nakamura’s expression of Eq. (3), it is not surprising they found remarkably good quantitative agreement between the two estimates.
4. Evidence from a numerical simulation
a. Lagrangian behaviors in the two coordinates
Although there are some inspiring works investigating particle motions in contour-based coordinates (e.g., Sallée et al. 2011; Shuckburgh and Haynes 2003), none of them relate the problems to the Y-coordinate random-work model as expressed by Eq. (4). Therefore, we will first compare particle behaviors in the two coordinates and show that particles in the Y coordinate follow the model of Eq. (4).
Figure 4 shows the snapshots of tracer distributions and particle positions collected within two regions R1 and R2 at the start of the fourth year, after the turbulence has fully developed. As an example, these two regions are selected with two considerations. First, the region should be small in both geographical and contour-based coordinates so that turbulent statistics within the regions remain approximately uniform. Second, the two regions cover the same tracer contour (say, the contour of 2.37 here), facilitating a clear demonstration of the along-contour variation of mixing regime in the present domain configuration. As a result, R1 is chosen to be close to the boundary jet (Fig. 2) and particles in R1 will be quickly injected into the frontal region where significant stirring processes, such as eddy shedding and eddy splitting, prevail. In contrast, R2 is chosen to be at relatively quiescent region where small-amplitude nonbreaking wave undulation prevails.
From Fig. 4 we can see that particles in R1 spread out quite well in x–y space after 250 days (Fig. 4f), which gives a false impression that particles are somewhat well mixed. These particles are only, if any, well mixed along tracer contours (Rhines and Young 1983) and their spreading in the cross-contour direction is expected to be extremely small as the tracer evolves quasi-adiabatically. One clear piece of evidence is that particles in R1 patch never drift to the upper half of the domain where tracer is shaded blue. For the R2 patch, as there is limited stirring effect, particles remain compact all the time and their spreading in both x–y and Y coordinates is also expected to be small.
To verify this, particle x and y positions are plotted in both coordinates in Fig. 5. Notice that all the position time series are shifted to facilitate the inspection of their dispersions. It is clear that in x–y coordinates (Figs. 5a,b), the R1 patch spreads gradually at the start and becomes two groups after ~50 days. When one group is advected into the frontal region at day 90, spreading becomes significant, especially in the x direction. Another group follows a similar way with a time lag of ~60 days. It is important to note that, as shown in both panels (Figs. 5a,b), particle motions in geographical space are highly correlated in time, as they might be trapped by looping currents or coherent vortices. Therefore, their velocity autocorrelation functions are expected to decay very slowly (to be verified later). For the R2 patch deployed in quiescent region (Figs. 5d,e), spreading is limited and all particles remain close to each other even after one model year. However, reversible undulation between 80 and 140 days is also evident, causing temporal correlated particle motions.
To demonstrate particle trajectories in the Y coordinate, one should first calculate the Y position of particles. There are two ways to do this. The first one is to use the precalculated relation to obtain the particle Y position given its tracer value q, and we call this the contour method. The second one follows the practice of Sallée et al. (2011) that obtains the Lagrangian changes of q between two consecutive times ( 1 day) and then estimate the cross-contour displacement as , where tracer gradient is evaluated in the Y coordinate. Finally, integrating from any constant start position (chosen to be 1000 km in Fig. 5) gives the position time series . In contrast to Sallée et al. (2011), here we choose to use the instant (rather than their climatological) local Y gradient to estimate the instantaneous . Practically, the gradient is an average between two consecutive times and hence
which facilitates the integration of in a trapezoidal fashion. We call this the gradient method.
The above two methods eventually provide very similar trajectory estimates in the Y coordinate (figure not shown). Therefore, only the results by the contour method are shown here as the precalculated is consistent with the transformed random-walk model of Eq. (4). Figures 5c and 5f show the particle trajectory in the Y coordinate for R1 and R2, respectively. Since the tracer is quasi-adiabatic, it is expected that the Y positions of all particles are close to horizontal lines. However, only those particles in R2 give satisfactory results (Fig. 5f). Particles in R1 (Fig. 5c), although parts of their tracks are close to straight lines, show considerable meridional cross-contour displacements (the spurious numerical diffusivity is causing the meridional motion in the Y coordinate). This is most significant after ~70 days when a group of particles in R1 drift into the frontal region where stirring effect is strong. In R2, particles also start to disperse around 90–130 days (Fig. 5f), corresponding to the stirring effect of a wave impulse (Figs. 5d,e). All these facts manifest that particle dispersion in the Y coordinate is subjected to small and enhanced by stirring.
Lagrangian velocity autocorrelation is also investigated (Fig. 6). In the x–y coordinates, a Eulerian time-mean flow is subtracted first before calculating the autocorrelation, following the standard practice (e.g., Bauer et al. 1998; Chiswell et al. 2007; Davis 1991; Poulain 2001; Qian et al. 2013). Velocity in the Y coordinate is obtained from through finite differencing, and the corresponding mean flow is assumed to be negligible. It is obvious that the Lagrangian velocities decorrelate very slowly in the x–y coordinates (Figs. 6a,b and 6d,e). Pronounced negative lobes are identified in R2 (Figs. 6d,e) and even after 50 days the autocorrelations do not tend to level off. This is associated with the back and forth meandering effect of nonbreaking waves (e.g., Klocker et al. 2012a). The Lagrangian decorrelation time scale, defined as the net area under the autocorrelation curves (e.g., Lumpkin et al. 2002), is roughly estimated as above 5 days for both regions. In contrast, the (mean) autocorrelations for both regions in the Y coordinate quickly level off within 5 days (Figs. 6c,f) and the first crossing-zero points are very close to the sampling interval of 1 day. This leads to a negligible Lagrangian decorrelation time in the Y coordinate. The very small negative lobs and oscillations could be attributed to 1) the mean drift in the Y coordinate is not removed here and particles are not averaged over the entire contour [only along-contour averaging recovers exactly Eq. (4)], and 2) practical errors induced by finite difference calculations. As the negative lobs is very small given that the time interval of the data is 1 day, the autocorrelations shown in Figs. 6c and 6f could be viewed as a delta function. Therefore, we conclude that particle motions in the Y coordinate follow the random-walk model and do decorrelate immediately in the contour-based coordinate, in significant contrast to those in geographic coordinates.
Figure 7 shows the different components of Lagrangian dispersion. The purpose here is to demonstrate different dispersion behaviors in the two coordinates rather than to interpret the mixing features, as particles are restricted within the model domain and shear dispersion induced by a time-mean flow is not removed. In Fig. 7, the x–y components of dispersion show significant spikes, indicating reversible stirring by waves, eddies or even basin-scale gyre circulation. In contrast, the Y-component dispersion in both regions remains close to zero until the enhancement of stirring at about day 80 for R1 and day 90 for R2. Another enhancement in R2 is also identified around day 240 (see Fig. 7b and Fig. 5f). In the rest of times, Y-component dispersion remains relatively flat, indicating quasi-conservative dynamics. In viewing the Lagrangian diffusivity, determined by the slope of dispersion in Fig. 7, one cannot expect any asymptotic behavior in x–y coordinates because particle spreading is restricted by boundaries and shear dispersion is present (e.g., Bauer et al. 1998). In the transformed Y coordinate, however, the growth rate of dispersion is already “asymptotic” at every time step. Note that the negative slope of Y dispersion in Fig. 7 does not prove the existence of negative diffusivity (or antidiffusion) because the Lagrangian diffusivity in the Y coordinate should be positively defined (the same as the effective diffusivity). This should be attributed to the inhomogeneity of diffusion [the drift correction term in Eq. (4)] as particles will sample different contours after a while from point release.
It is worth noting that, in the present model configuration, neither the effective diffusivity nor the traditional particle diffusivity is useful to access the local mixing, which is often needed for Eulerian numerical models. It is already known that effective diffusivity only represents the along-contour mean diffusion (Nakamura 2008), and flow regimes demonstrated here could vary substantially along the contour (comparing R1 and R2). In the Lagrangian framework, particles generally spread out and sample different flow regimes as time goes on and it is hard to interpret the averaged results (e.g., Fig. 7) at large time lags. In this case, one can use the instantaneous Lagrangian dispersion rate as the asymptotic diffusivity based on the fact that particle decorrelates very quickly in the Y coordinate (Figs. 6c,f). This also provides an opportunity to demonstrate the full temporal evolution of Lagrangian diffusivity, which is generally impossible in geographic coordinates if one has to pursue the asymptotic one. Such an instantaneous Lagrangian dispersion rate in the Y coordinate can also be used to verify its equivalence to the effective diffusivity, as given in the following subsection.
b. Evidence for the equivalence of the two diffusivities
Lagrangian diffusivity will be estimated using Eq. (6), the half growth rate of particle dispersion. Nakamura’s effective diffusivity is calculated as
where the operator is understood as an along-contour average and will be evaluated using a finite-difference approximation (Nakamura 2008). Here we choose as the minimum possible length (it is not actually minimum as ) because it renders as a meridional diffusivity in the Y coordinate and the unification of the two diffusivities does not depend on this choice. A 4-yr-average value of 3 m2 s−1 (Fig. 3) will be used in Eq. (7).
As a first step to the reconciliation of the two diffusivities, the first moment of particle displacement (i.e., mean drift) will be presented. This is equivalent to the drift correction term in Eq. (4). To show this, particles along a tracer contour q at an instant time t are collected as an ensemble patch. The displacement of each particle is then estimated as , where is one day and particle position Y will be estimated using both the contour method and gradient method described in section 4a. Through ensemble averaging and repeating over all the contours and time steps, the mean drift , or , can be obtained. As a reference, the mean drift is also estimated using Eq. (5) with 1 day.
All these estimates are shown in Fig. 8. From Fig. 8a one can see that elevated mean drift emerges after 120 days and spreads quickly toward the boundaries. After about one year, two obvious mean drift bands, one around 1800 km with negative values and one around 300 km with positive values, become robust and move slowly toward the boundaries. Two separating zero lines corresponds to the maxima of effective diffusivity (Fig. 9a), and the mean drift bands converge into two regions of strong mixing (e.g., Davis 1991). There are also alternating thin drift bands between the two. Among these it is worth noting that around the midbasin line ( 1100 km), there is a dipole-like concentrated drift bands, indicating a convergence at the midbasin line. All these gross features are well reproduced by the mean drift of Lagrangian particles, as demonstrated in Figs. 8b and 8c using different methods.
Now we move on to the second moment of particle displacement, that is, Lagrangian dispersion diffusivity. The reference is the Nakamura’s effective diffusivity (Fig. 9a) obtained using Eq. (7). The evolution of tracer contours, as superimposed in Fig. 9a, are relatively straight and only disperse slowly toward the boundaries, indicating quasi-adiabatic dynamics. The tracer gradient in the Y coordinate is weak where the effective diffusivity is high while tracer front (strong gradient indicated by dense contours) occurs in regions of low diffusivity. This is because tracer contours dispersed faster in regions of high effective diffusivity and remain compact where diffusivity is weak. It is interesting to note that there is a slight densification of contours between 1.9 and 2.1 during 240–480 days. This slight sharpening of tracer gradient cannot be viewed as antidiffusion or upgradient transport. In the Y coordinate, this is just caused by inhomogeneity of effective diffusivity [i.e., the drift correction term in Eq. (4)] and eventually all tracer gradient will be irreversibly removed.
In addition to Nakamura’s standard method, we also devise an intermediate estimate using Lagrangian data. This method is used to check whether Lagrangian sampling (over 1000 particles per contour) is sufficient to reproduce a similar pattern as in Fig. 9a. For this, we write Eq. (7) as
Note that is a 2D field and is a 1D field. The along contour average operator in Eq. (7) is exactly the same as the ensemble average operator taken over all the Lagrangian particles along a single contour q. Gradient operators and are defined respectively in each coordinate (). The second equality holds because of (i.e., does not change along the contour q). The third equality holds because of . We first calculate the gradients and interpolate them linearly to particle positions in the two coordinates to obtain the gradient ratio . The gradient ratio is then squared, ensemble averaged over all the particles along each contour and finally multiplied by 3 m2 s−1, to obtain another estimate of effective diffusivity. This method is fundamentally identical to Nakamura’s method [i.e., Eq. (7)] and the differences lie in that the gradient ratio is sampled by Lagrangian particles or by Eulerian grid points. Figure 9b shows this estimate, which resembles the reference shown in Fig. 9a except that the overall magnitude is somewhat smaller. Such a fact has confirmed that more than 1000 particles per contour are sufficient to provide robust estimates.
To demonstrate the full temporal evolution, the Lagrangian dispersion diffusivity, defined by the left-hand-side of Eq. (6), will be practically discretized as
where and 1 day. This only requires particle Y positions at two neighboring times. As one has to recollect particles along a contour at every time instant t, the collected ensemble particles can be treated as newly “point released” in the Y coordinate at time t; they are actually along-contour collected in x–y coordinates. Besides, the commonly used pseudotrack method (e.g., Griesel et al. 2010; Klocker et al. 2012b; Swenson and Niiler 1996) is not applicable here as it involves a time average. As a result, there are ~1000 independent particles at every grid point in the Y coordinate. It is worth noting that as the pseudotrack method cannot be used, the number of particles required for a robust estimate should be much larger than that of the practice in the real ocean. Tracking synthetic particles in numerical models is a good way to remedy this (e.g., Griesel et al. 2010; Klocker et al. 2012b; LaCasce et al. 2014).
The Lagrangian diffusivity estimates shown in Figs. 9c and 9d, respectively use the gradient method and the contour method to obtain the Y position. Previous estimates using the pseudotrack method only provided a climatological diffusivity. Now the full temporal evolution of the Lagrangian diffusivity in addition to its spatial distributions (in the Y coordinate) is available under the present analysis framework. This can be achieved only when instantaneous dispersion rate is already asymptotic. Comparing to Fig. 9a, it can be seen that there is a remarkably good quantitative agreement between the Lagrangian diffusivity and Nakamura’s effective diffusivity, not only in their spatial distributions but also in their temporal evolutions. The raw Lagrangian diffusivity estimates are a little noisy because the finite-time differencing introduces gridscale noise. A slight 5-day temporal running mean is used to suppress this gridscale noise and provides better results (Figs. 9c,d). A careful comparison reveals that the elevated diffusivity at the midbasin line ( 1100 km), which is obvious in Figs. 9a and 9b, is somewhat underestimated by the Lagrangian diffusivity. Although the gradient method yields a slightly weak diffusivity (Fig. 9c) than the contour method does (Fig. 9d), both estimates reproduce the overall features, as well as specific details such as the occurrence time and location of diffusivity maxima/minima. All in all, the Lagrangian diffusivity in the Y coordinate does reproduce many of the features in Fig. 9a. All the evidence from Figs. 8 and 9 verifies that Lagrangian particle motions in the transformed Y coordinate do satisfy the random-walk model of Eq. (4). More importantly, Lagrangian diffusivity is exactly effective diffusivity only in the Y coordinate, rather than in the geographic x–y coordinates.
5. Conclusions and discussion
Although there are many ways to access the temporal and spatial variations of lateral eddy diffusivity, different definitions for the diffusivity has hampered a consistent diffusivity estimate and a physically based interpretation. This leads to the present reconciliation of two diffusivities: the Lagrangian diffusivity and the effective diffusivity. Through using the MITgcm model, a dynamical flow is generated to drive the evolutions of a quasi-conservative tracer q and a large number of Lagrangian particles. The tracer contours are then used to construct a new spatial Y coordinate, and in this contour-based Y coordinate, theoretical and practical discrepancies between the Lagrangian diffusivity and effective diffusivity are well reconciled. In the Y coordinate, it is shown that the Lagrangian particles obey the transformed zeroth-order stochastic model of Eq. (4), whose random impulse is exactly regulated by the effective diffusivity. The Fokker–Planck equation of the random-walk model is the diffusion-only equation already proposed by Nakamura (1996) in contour-based coordinates. Estimates of the first and second moments of particle displacements, as well as the transport experiments done by Shuckburgh and Haynes (2003, see their Fig. 12 and relevant discussion), all support that Eq. (4) does control the particle motion in the transformed Y coordinate.
Previous issues unaddressed by theoretical treatment of Lagrangian and effective diffusivities include 1) the dispersion-based diffusivity definition does not have an explicit connection to the small-scale diffusion (e.g., Vallis 2006); 2) the time used for assessing the Lagrangian diffusivity must be past the decorrelation time scale, which eventually hampers localized estimates at a time instant; 3) previous reconciliation may depend somewhat on the zonal periodicity of the domain and effective diffusivity is not strictly guaranteed to be in the meridional direction even when zonal periodicity is ensured (e.g., Wolfram and Ringler 2017); and 4) the dispersion-based diffusivity is a tensor while the effective diffusivity is a scalar. Now these issues can be well addressed as 1) dispersion-based diffusivity does depend on small-scale diffusion. In the absence of small-scale diffusivity (e.g., no molecular diffusivity), particles do not disperse at all in the Y coordinate and Lagrangian diffusivity becomes zero (Fig. 1). In this sense, Lagrangian diffusivity can be directly linked to the irreversible mixing; 2) the instantaneous particle spreading rate can be readily treated as the asymptotic Lagrangian diffusivity, as Eq. (4) is a zeroth-order stochastic model (not higher-order ones). Skipping the decorrelation time scale is naturally satisfied and estimates are local (in the Y coordinate) and instantaneous; 3) the Y coordinate is more general and does not require periodicities for the domain because the “zonal” mean in the Y coordinate is exactly the along-contour mean; and 4) the Lagrangian diffusivity tensor reduces naturally to a scalar in the Y coordinate, corresponding to the effective diffusivity that measures only cross-Y-surface (irreversible) mixing.
As the above issues can now be well addressed based on Eq. (4), we conclude that the two diffusivities have been unified, and can only be exactly unified in contour-based coordinates. One may ask the question that can they be unified in Eulerian (geographic) x–y coordinates; or put it another way, whether it is possible to reconciling Eulerian diffusivity and effective diffusivity. The answer is no, unfortunately, or at least not straightforward to our current knowledge. This is because the core benefit of effective diffusivity concept is to separate the effect of irreversible mixing from that of adiabatic stirring. To achieve this, one has to abandon the Eulerian positional information and measure fluid motion relative to conservative tracer contour, that is, the Y-coordinate surfaces, rather than Eulerian x–y surfaces (Nakamura 1998). As a result, one can sense the irreversible mixing by observing the cross-contour motion (but also lose the ability to obtain the along-contour motion, that is, cannot quantify the anisotropy of mixing). If, taking a step back, there exists a way of reconciling Eulerian diffusivity and effective diffusivity, it should also be built on the physical basis that adiabatic stirring should be isolated from irreversible mixing (Table 1). Studies (e.g., Eden et al. 2007, 2009; McDougall and McIntosh 2001; Medvedev and Greatbatch 2004; Nakamura 2001; Osborn and Cox 1972) have already shed some lights along this path.
Traditional calculation of the Lagrangian diffusivity generally involves a scale-separation process () and uses the eddy part to infer Lagrangian diffusivity. As advective stirring, or shear dispersion, can be caused by motion at any scales (Winters and D’Asaro 1996), shear-augmented dispersion by still exists even when is removed. Therefore, contributions from shear dispersion will be inevitable. The magnitude of Lagrangian diffusivity is then determined by shear dispersion of (small-scale diffusion only makes a tiny contribution). In the limit that the separation scale shrinks to the viscous scale or the grid size of a numerical model, Lagrangian diffusivity becomes close to the small-scale diffusivity and hence mostly measures irreversible mixing. That is, the more we put into , the more the shear-augmented effect will likely be, and the larger estimates one will get. In the perspective of measuring irreversible mixing, an adiabatic shear contribution should be removed from Lagrangian dispersion. The concept of effective diffusivity inherently uses the model gridscale to accomplish scale separation. As a result, it avoids any explicit, which is also somewhat arbitrary, decomposition of a flow into a mean and associated eddy part (Nakamura 1998). Lagrangian diffusivity, if it incorporates an Eulerian bin to accomplish scale separation, will become a mixed Eulerian–Lagrangian diagnostic (LaCasce 2008) and lose its Lagrangian nature to isolate the adiabatic stirring (shear) from irreversible mixing.
The present method, calculating dispersion using displacement across instantaneous contours, can also be applied to real oceans. One may combine tracer contour information (sea surface height or temperature/salinity) and surface drifters (Lumpkin and Pazos 2007) to infer Lagrangian dispersion in contour coordinates, as already practiced by some inspiring works (e.g., Roach et al. 2016; Sallée et al. 2011). However, at the sea surface, diabatic processes (e.g., heat and freshwater fluxes) other than diffusion will certainly contaminate the diffusivity estimates. One solution, similar to previous studies (e.g., Abernathey and Marshall 2013; Klocker and Abernathey 2014; Marshall et al. 2006; Shuckburgh et al. 2009b), is to drive a passive tracer whose evolution is subject only to diffusion with satellite-based surface currents. In the deep ocean, potential temperature and salinity may be good choices for seeding Lagrangian particles as the ocean interior is almost adiabatic (e.g., Ledwell et al. 2011; Wu et al. 2011). Direct estimates of in situ irreversible mixing during field campaigns where tracer patches are released in combined with some Lagrangian observations (e.g., Banyte et al. 2013; Boland et al. 2015; Tulloch et al. 2014), are also possible (e.g., Lee et al. 2009) given that Lagrangian floats can sample the tracer simultaneously. This will be practiced and reported in a future study.
This work was jointly supported by National Natural Science Foundation of China (Grants 41521005, 41676016), the CAS/SAFEA International Partnership Program for Creative Research Teams, Science and Technology Program of Guangzhou (Grant 201607020043), and Science and Technology Planning Project of Guangdong Province (Grant 20150217). The authors gratefully acknowledge the two anonymous reviewers for their constructive comments, as well as the use of the HPCC at the South China Sea Institute of Oceanology, Chinese Academy of Sciences.