Reconstruction of Diffusion Coefficients and Power Exponents from Single Lagrangian Trajectories

Abstract

Using discrete wavelets, a novel technique is developed to estimate turbulent diffusion coefficients and power exponents from single Lagrangian particle trajectories. The technique differs from the classical approach (Davis (1991)’s technique) because averaging over a statistical ensemble of the mean square displacement (<X²>) is replaced by averaging along a single Lagrangian trajectory X(t) = {X(t), Y(t)}. Metzler et al. (2014) have demonstrated that for an ergodic (for example, normal diffusion) flow, the mean square displacement is <X²> = $\underset{T\to \infty }{\lim }{\tau }_X^2\left(\mathrm{T},\mathrm{s}\right)$, where ${\tau }_X^2$ (T, s) = 1/(T − s) ${{\displaystyle \mathsf{\int }}}_0^{T-s}(X\left(t+\Delta t\right)$ − X(t))² dt, T and s are observational and lag times but for weak non-ergodic (such as super-diffusion and sub-diffusion) flows <X²> = $\underset{T\to \infty }{\lim }\ll {\tau }_X^2\left(\mathrm{T},\mathrm{s}\right)\gg$, where $\ll \dots \gg$ is some additional averaging. Numerical calculations for surface drifters in the Black Sea and isobaric RAFOS floats deployed at mid depths in the California Current system demonstrated that the reconstructed diffusion coefficients were smaller than those calculated by Davis (1991)’s technique. This difference is caused by the choice of the Lagrangian mean. The technique proposed here is applied to the analysis of Lagrangian motions in the Black Sea (horizontal diffusion coefficients varied from 10⁵ to 10⁶ cm²/s) and for the sub-diffusion of two RAFOS floats in the California Current system where power exponents varied from 0.65 to 0.72. RAFOS float motions were found to be strongly non-ergodic and non-Gaussian.

1. Introduction

Ivanov and Chu [1] proposed an original method for the analysis of drifters/floats deployed in the ocean as single Lagrangian particles. Following that paper, the discrete wavelet decomposition (the 10th order Daubechies wavelet is used in the present paper) [2,3] was applied to observed drifter/float trajectories Xⁱ = (Xⁱ, Yⁱ), i = 1,…I (where I is the total number of drifters/floats), and found that the real drifter/float motions can be divided into two parts: the mean Lagrangian (regular) motions

X_mean (X_o, Y_o, t) ≅ X₁ + X₂ +…+ X_P,

(1)

and irregular components, which may be random as well as deterministic chaos,

X_irreg (X_o, Y_o, t) ≅ X_{P + 1} + X_{P + 2}+…+ X_L,

(2)

| Xobs − Xmean-Xirreg | ≤ ϵ(L),

(3)

ϵ(L)→ϵ_o, if L→∞,

(4)

where ϵ_o is an unlimited small value.

Here it is assumed that the irregular component (2) may be described as the diffusion of drifters/floats, which is determined through the tensor of diffusion coefficients (K_ij). Instead of averaging over a statistical ensemble, averaging along a Lagrangian trajectory is introduced and the ergodic hypothesis to used estimate the diffusion coefficients. This allows the single Lagrangian trajectories to be analyzed and turbulent diffusion coefficients estimated (or other transport coefficients if needed) along the single Lagrangian trajectory. A number of techniques have been applied to decompose Lagrangian trajectories into their mean and diffusive components, for example, by binning velocities and then fitting them with cubic splines [4,5] or by applying a Gauss-Markov estimator [6,7]. It was pointed out [8] that some of these techniques are based on making assumptions about the correlation between Lagrangian and Eulerian representations. Others [4] used datasets too small to allow for checking their convergence. Besides, as discussed in [8,9], there is an ambiguity in this decomposition, with different techniques giving different results.

Techniques to estimate turbulent diffusion coefficients from Lagrangian trajectories have been developed by [5,10,11,12,13,14,15,16] among others, see [1] for more details. In this paper a novel technique for estimation of turbulent diffusion coefficients from a single Lagrangian trajectory has been developed and will subsequently be referred to as the Single Trajectory Method or “STM.” The STM will be applied to the analysis of 13 drifters deployed in the Black Sea in 2002–2003 [17] and two floats deployed in the California Current System [18]. In the latter case, floats were chosen which demonstrated abnormal eddy diffusion (sub-diffusion) caused by the trapping of floats by mesoscale eddies.

Since the drifters and floats can be described as single particles, usual statistical ensembles are impossible to construct. Therefore, the theoretical ideas from [19,20] where one-dimensional examples are discussed as well as the theoretical paper of [21] who analyzed a multi-dimensional model with statistically variable diffusion coefficients. The general idea from [19,20] is to replace the averaging over a statistical ensemble by averaging along a single trajectory of a drifter/float. This idea was developed for both pure ergodic and weak non-ergodic flows. This approach can be used in two cases: the analysis of a single drifter/float trajectory (examples are given in this paper), and for analysis of the ensemble of Lagrangian trajectories to reconstruct each as a single trajectory. Once the trajectories have been reconstructed, a spatial interpolation procedure is applied to re-calculate the results into nodes of a regular grid (see, for example, [22]). Note that the procedure does not require the turbulent field of X_irreg to be differentiable.

There are important differences between Ivanov and Chu’s (2019) approach and the one introduced in the present paper. A novel, highly-effective method to estimate the turbulent diffusion coefficients using the time-averaged mean square displacements has been developed, along with a filter (Equation (5)) for selecting the mean Lagrangian component and reducing the degree of wavelets to improve the resolution of the results. An original technique for estimating errors for diffusion coefficients and power exponents has been introduced, as described in Appendix B. The methods developed in the present paper allowed us to formulate a special approach for the analysis of individual Lagrangian trajectories, which can be applied to analyze regular as well as irregular and quasi-random processes. The present paper uses new data to calculate diffusion coefficients (drifter trajectories in the Black Sea) and power exponents (RAFOS floats deployed in the California Current system); anisotropy of turbulent diffusion has been estimated for the Black Sea.

The general methods for decomposition of a trajectory into the mean Lagrangian flows and turbulent diffusion will be discussed in Section 2. A simple parameterization for turbulent diffusion coefficients is described in Section 3. Section 4 contains discussions of the theoretical methods for the analysis of single Lagrangian particles as the so-called time averaged characteristics. Lagrangian data are analyzed in Section 5. Different numerical methods for the estimates of turbulent diffusion coefficients are discussed in Section 6. Turbulent diffusion was estimated from thirteen Lagrangian drifter trajectories in the Black Sea and two RAFOS floats in the California Current system. In the latter case, degrees of non-ergodicity and non-Gaussivity for the floats are discussed and power exponents for sub-diffusion are calculated. Details are given in Section 7. Section 8 contains conclusions. Appendix A and Appendix B discuss the moments of time averaged characteristics of Lagrangian particles and how errors for determination of turbulent diffusion coefficients and power exponents can be found.

2. Decomposition into Mean Lagrangian and Turbulent Trajectories

A Lagrangian trajectory X = (X, Y) is decomposed using a Daubechies wavelet of the L-order to divide the observed Lagrangian trajectory X into L sub-trajectories X₁ (τ₁),…, X_L (τ_L) where τ₁ > τ₂…> τ_L are some specified time scales and X $\cong$ X₁ + X₂ +…+ X_L. The last approximation will be close to the exact one when the data quality is high (measurement errors are less than 0.5%). First, the mean Lagrangian flow is determined as X_mean ≅ X₁ + X₂ +…+ X_P and specification of P is discussed assuming a regular character of the mean Lagrangian flow (if necessary).

This is done from analysis of the phase spectrum introduced by the method of [23] for each sub-trajectory and provides a ratio of the imaginary and real parts of the Fourier transform of this trajectory. A regular trajectory has minimal possible degree of filling of this phase space [1]. A boundary between the mean Lagrangian and turbulent flows is given a priori. It is necessary to estimate P.

If noise contributions to the mean Lagrangian flow are quite strong, the noise needs to be filtered out. A recursion filter [24] is applied to reduce the noise contribution. It also determines the nature of observed irregular behavior. Both enable one to construct a new trajectory in which the effective noise contribution is lower.

Here we describe the recursion filter developed by [24]. If we have an observation sample X = {X₁(t),…, X_{k + 1}(t)} measured at points t = t₁,…, t_{k + 1} then

$$\mathrm{Yr}\left({\mathrm{t}}_{\mathrm{k} + 1}\right) = \frac{{\sum }_{\mathrm{n}=1}^{\mathrm{N}}\mathrm{s}(\mathsf{\epsilon }{¯}^1|X\left({\mathrm{t}}_{\mathrm{n}}\right)-Y\left({\mathrm{t}}_{\mathrm{k}}\right)| {\mathrm{X}}_{\mathrm{r}}\left({\mathrm{t}}_{\mathrm{n}+1}\right)}{{\sum }_{\mathrm{n}=1}^{\mathrm{N}}\mathrm{s}(\mathsf{\epsilon }¯¹\left|X\left({\mathrm{t}}_{\mathrm{n}}\right)-Y\left({\mathrm{t}}_{\mathrm{k}}\right)\right|)}, \mathrm{r} = 1,2$$

(5)

where ε << 1 is a small parameter, $s\left(u\right)$ is a weighting function which tends to zero when $u\to \infty , \mathrm{and} N \mathrm{is} \mathrm{the} \mathrm{number} \mathrm{of} \mathrm{observations}$.

It was demonstrated in [25] that the best option here is when $s=\exp \left(-u^2\right)$ and ε~0.01. The variance of a random signal is proportional to 1/Nε and it tends to zero when Nε$\to \infty$ for ε$\to 0$ and N$\to \infty .$ The deterministic part of the signal, no matter how close the final is to deterministic chaos, has an additional term proportional to ε², which tends to zero when ε $\to 0$. The proposed method does not require the noise to be small. Note also that deterministic filtering models (5) can be applied to regular as well as irregular and quasi-random processes.

Our numerical experiments demonstrated that the mean Lagrangian flow along a single trajectory is a complex function of time, and the non-dimensional parameter E_mean/E_turb determines the value of turbulent diffusion [9] where E_mean and E_turb are the energy of the mean Lagrangian and turbulent flows averaged along the Lagrangian trajectory, respectively.

Then, we determine the residual δX = X − X_mean and analyze the behavior of δX along the Lagrangian trajectory assuming that irregular motions can be described by the diffusion tensor coefficients. In principle, other methods for description of diffusion motions exist (see, for example, [26]) but are not described here.

3. Simple Diffusion Parameterization

An approach [27] was applied to estimate the two diagonal diffusion coefficients (K₁₁ and K₂₂) proportional to the mixing length (λ) and the intensity of turbulent velocity (λ ε)^1/2 [28], where ε is the rate of energy dissipation. Note that other methods (for example, [26]) can be used as well.

Symmetrical coefficients of tensor ${\mathrm{K}}_{12}^s$ and ${\mathrm{K}}_{12}^a$ are obtained from the dimensionality ideas and definitions which can be found in [29] (see, [1,27]):

$${\mathrm{K}}_{11} = \mathrm{c} {\mathsf{\sigma }}_1 {\mathsf{\delta }}_1, {\mathrm{K}}_{22} = \mathrm{c} {\mathsf{\sigma }}_2{\mathsf{\delta }}_2,{ \mathrm{K}}_{12}^s = \mathrm{c} \left({\mathsf{\sigma }}_1 {\mathsf{\delta }}_2 + {\mathsf{\sigma }}_2 {\mathsf{\delta }}_1\right),$$

(6)

where c = 0.1 [30], ${\mathsf{\delta }}_1=\sqrt{1/\left(N-1\right)\sum _{\mathrm{n}=1}^{\mathrm{N}}{\mathsf{\delta }\mathrm{X}}_{\mathrm{n}}^2}$, ${\mathsf{\sigma }}_1= \sqrt{1/\left(N-1\right)\sum _{\mathrm{n}=1}^{\mathrm{N}}{\mathsf{\delta }\mathrm{U}}_{\mathrm{n}}^2}$, ${\mathsf{\delta }}_2=1/\left(\mathrm{N}-1\right)\sqrt{\sum _{\mathrm{n}=1}^{\mathrm{N}}{\mathsf{\delta }\mathrm{Y}}_{\mathrm{n}}^2}$, ${\mathsf{\sigma }}_2=\sqrt{1/\left(N-1\right)\sum _{\mathrm{n}=1}^{\mathrm{N}}{\mathsf{\delta }\mathrm{V}}_{\mathrm{n}}^2}$, {δU_n, δV_n} is the residual velocity along the single Lagrangian trajectory, $K^s$ is the symmetric tensor of K.

Numerical experiments with the Black Sea drifters, RAFOS floats and drifters deployed in the California Current System [1] have demonstrated that a turbulent trajectory is differentiable if we exclude rapid motions which contribute very little to the total Lagrangian trajectory.

To use equations for determination of diffusion coefficients, averaging is introduced because often the coefficients do not converge for t→∞ or converge asymptotically ($\left|\right|\mathrm{K}\left|\right|$ is limited by the upper K_upper and lower K_low boundaries and changes with time within [K_upper, K_low]). Assuming that the Lagrangian trajectory is quite long, averaging over an ensemble is replaced by averaging along this trajectory.

4. Methods for Calculation of Diffusion Coefficients

For a single Lagrangian trajectory (or when the number of trajectories is not high), diffusion coefficients K can be deduced in another way (see, for example, [19,20]) than in classical works (for example, [31]). This new approach uses the definition of ergodicity and replaces the averaging over a statistical ensemble by the averaging along the Lagrangian trajectory.

Following [19,32], the time averaged mean squared displacements τ_X² (T, s) are introduced for the residual of a Lagrangian trajectory δX = X − X_mean and other moments of higher order for ${\tau }_X^2$ (T, s) as

$${\tau }_X^{\mathrm{q}}\left(\mathrm{T}, \mathrm{s}\right) = 1/\left(\mathrm{T} - \mathrm{s}\right){{\displaystyle \mathsf{\int }}}_0^{\mathrm{T}-\mathrm{s}}{\left(\mathsf{\delta }X\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }X\left(\mathrm{t}\right)\right)}^{\mathrm{q}} \mathrm{dt},$$

(7)

If q = 2 the following equations are correct:

$${\tau }_{\mathrm{X}}^2\left(\mathrm{T}, \mathrm{s}\right) \approx 1/\left(\mathrm{T} - \mathrm{s}\right){{\displaystyle \mathsf{\int }}}_0^{\mathrm{T}-\mathrm{s}}\left(\mathsf{\delta }\mathrm{X}\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }\mathrm{X}\left(\mathrm{t}\right)\right) \left(\mathsf{\delta }\mathrm{X}\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }\mathrm{X}\left(\mathrm{t}\right)\right) \mathrm{dt},$$

(8)

$${\tau }_{\mathrm{Y}}^2\left(\mathrm{T}, \mathrm{s}\right) \approx 1/\left(\mathrm{T} - \mathrm{s}\right){{\displaystyle \mathsf{\int }}}_0^{\mathrm{T}-\mathrm{s}}\left(\mathsf{\delta }\mathrm{Y}\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }\mathrm{Y}\left(\mathrm{t}\right)\right) \left(\mathsf{\delta }\mathrm{Y}\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }\mathrm{Y}\left(\mathrm{t}\right)\right) \mathrm{dt},$$

(9)

$${\tau }_{\mathrm{XY}}^2\left(\mathrm{T}, \mathrm{s}\right) \approx 1/\left(\mathrm{T} - \mathrm{s}\right){{\displaystyle \mathsf{\int }}}_0^{\mathrm{T}-\mathrm{s}}\left(\mathsf{\delta }\mathrm{X}\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }\mathrm{X}\left(\mathrm{t}\right)\right) \left(\mathsf{\delta }\mathrm{Y}\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }\mathrm{Y}\left(\mathrm{t}\right)\right) \mathrm{dt},$$

(10)

where T and s are the observational period and lag time, respectively.

Since the residual Lagrangian trajectory describes only a pure diffusion process without mean Lagrangian drift, a simple diffusion model [21,33] can be used to describe such behavior

dδX_i/dt = (2)^1/2D_ij(t) f_j(t), <f_i> = 0, <f_i(t)f_j(t’)> = δ_ij δ(t − t’), i, j = 1,2

(11)

where δ is the delta function, δ_ij = 0(1) if i≠j (i = j), to determine

$${\tau }_X^2\left(\mathrm{T},\mathrm{s}\right) = 1/\left(\mathrm{T} - \mathrm{s}\right) {{\displaystyle \mathsf{\int }}}_0^{\mathrm{T}-\mathrm{s}}{\left(\mathsf{\delta }X\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }X\left(\mathrm{t}\right)\right)}^2 \mathrm{dt} = 2 \mathrm{Sp} \left(K\right) \mathrm{s},$$

(12)

where K = DD^T is the tensor of turbulent diffusion coefficients.

Equation (15) was analytically solved by [21]. Here integrals in the solutions were determined when diffusion coefficients were neither stochastic nor developing from the time t.

It was demonstrated [19,20] that, for ergodic flows (for example, the normal diffusion)

$$\langle X^2\rangle = \underset{T\to \infty }{\lim }{\tau }_X^2\left(\mathrm{T}, \mathrm{s}\right),$$

(13)

i.e., the mean obtained for a statistical ensemble is equal to the mean along a single trajectory.

For weak non-ergodic (super-diffusion and sub-diffusion) flows, this equation should be replaced by

$$\langle X^2\rangle = \underset{T\to \infty }{\lim }\ll {\tau }_X^2\left(\mathrm{T}, \mathrm{s}\right)\gg ,$$

(14)

where additional averaging over an L-dimensional statistical ensemble is introduced as

$$\begin{array}{cc}\hfill \ll {\tau }_X^2\left(\mathrm{T}, \mathrm{s}\right)\gg = 1/L{\sum }_{\mathrm{l}=1}^{\mathrm{L}}{\left\{\mathsf{\tau }{\mathrm{x}}^2\left(\mathrm{T}, \mathrm{s}\right)\right\}}_{\mathrm{l}} =& 1/\left(\mathrm{T} - \mathrm{s}\right){{\displaystyle \mathsf{\int }}}_0^{\mathrm{T}-\mathrm{s}}<{\left(\mathsf{\delta }\mathbf{X}\left(\mathrm{t}+\mathrm{s}\right)-\mathsf{\delta }\mathbf{X}\left(\mathrm{t}\right)\right)}^2> \mathrm{dt},\hfill \\ \hfill \mathrm{l} =& 1,\dots ,\mathrm{L};\hfill \end{array}$$

(15)

Calculations under the conditions K₁₂ = K₂₁ and $\parallel \mathrm{K}\parallel \ge 0$(assuming that antisymmetric diffusion can be neglected because its coefficients are very small in comparison to K₁₁ and K₂₂ and the drifters/floats are moving in the deep parts of the ocean/sea) give the simplest equations for determination of elements of the tensor of turbulent diffusion coefficients as

$${\mathrm{K}}_{11} \approx 1/2 \mathrm{s} {\tau }_{\mathrm{X}}^2\left(\mathrm{T}, \mathrm{s}\right),$$

(16)

$${\mathrm{K}}_{22} \approx 1/2 \mathrm{s} {\tau }_{\mathrm{Y}}^2\left(\mathrm{T}, \mathrm{s}\right),$$

(17)

$${\mathrm{K}}_{12} \approx 1/2 \mathrm{s} {\tau }_{\mathrm{XY}}^2\left(\mathrm{T}, \mathrm{s}\right).$$

(18)

In principle there is another way to estimate K₁₂ from the analysis of the fourth moment for ${\tau }_X$ (T, s) if the Lagrangian trajectory is long enough:

$$\langle \langle {\tau }_X^4\left(\mathrm{T},\mathrm{s}\right)\rangle \rangle = 4 \times {\mathrm{s}}^2 \times \left[\mathrm{Sp}\left(K\right) \times \mathrm{tr}\left(K\right) + 2 \times \mathrm{Sp}\left(K \times K\right)\right].$$

(19)

Simple calculations give

$${\mathrm{K}}_{12} = {\{1/2 \left[\langle \langle {\tau }_X^4\left(\mathrm{T},\mathrm{s}\right)\rangle \rangle - 4 \times {\mathrm{s}}^2 \times \mathrm{Sp}\left(K\right) \times \mathrm{Sp} \left(K\right)\right]/\left(8 \times {\mathrm{s}}^2\right) - {\left({\mathrm{K}}_{11}\right)}^2 - {\left({\mathrm{K}}_{22}\right)}^2\}}^{1/2}.$$

(20)

Note that Equations (19) and (20) are correct if D does not depend on time.

To account for the weak non-ergodicity of flows, an additional averaging <<…>> is introduced into Equation (7) and <<${\tau }_X^2$ (T, s)>> is re-calculated. Here simple models were used to obtain equations for sub-diffusion dynamics which describes the trapping of particles by mesoscale eddies.

For example, <<${\tau }_X^2$ (T, s)>> is calculated in a multi-dimensional case assuming that X is described by a continuous time random walk model. This was obtained [34] using the analytical solution for a multi-dimensional fractal equation and found that

$$\langle \langle {\tau }_X^2\left(\mathrm{T},\mathrm{s}\right)\rangle \rangle ~{\mathrm{s}}^{1-\mathsf{\alpha }}, \mathrm{for} \mathrm{t} >> s (\mathrm{sub}-\mathrm{diffusion} \mathrm{behavior}, 0 < \alpha \le 1),$$

(21)

$$\langle \langle {\tau }_X^2\left(\mathrm{T}, \mathrm{s}\right)\rangle \rangle ~\mathrm{s}, \mathrm{for} \mathrm{t} << \mathrm{s}, \left(\mathrm{normal} \mathrm{diffusion}, \alpha = 0\right),$$

(22)

which agrees with results of [19,20] to within a constant in Equations (21) and (22). These results differ from conclusions [18] which showed that ballistic diffusion ($\alpha$ = 2) should be observed before sub-diffusion. These conclusions are checked below using two RAFOS floats as an example.

Note that if Equations (21) and (22) describe the case of sub-diffusion, for large T the long-time average will provide the same information as the ensemble average. The general theory for ergodic and weak non-ergodic flows was developed [19] and applied the theory for sub-normal diffusion processes [20].

Where do ergodicity and Gaussivity not hold? For ergodicity, this can be determined using the variance of the amplitude fluctuations of ${\tau }_X^2$ (T, s) [19]

$$\mathrm{EB} = \underset{t\to \infty }{\lim } [\langle {\xi }^2\rangle - {\langle \xi \rangle }^2], \mathsf{\xi } = {\tau }_X^2\left(\mathrm{T}, \mathrm{s}\right)/\langle \langle {\tau }_X^2\left(\mathrm{T}, \mathrm{s}\right)\rangle \rangle ,$$

(23)

where EB is the degree of ergodicity and varies from EB = 1 for $\alpha =1$ to EB = 0 for $\alpha$ = 0 (Gaussian process). The non-Gaussivity parameter is determined as suggested by [20]

$$\mathrm{H} = 1/2 \langle \langle {\tau }_X^4\left(\mathrm{T},\mathrm{s}\right)\rangle \rangle /{\langle \langle {\tau }_X^2\left(\mathrm{T},\mathrm{s}\right)\rangle \rangle }^2 - 1,$$

(24)

and becomes non-zero if the distribution of the displacement (the van Hove correlation function) is not Gaussian.

According to a classical definition of ergodicity used in theoretical physics, the mean over time should converge to the mean over a statistical ensemble for an ergodic process when T$\to \infty$. However, we also know that any pure diffusion process characterized by a zero velocity and non-zero diffusion coefficients is an ergodic process. Therefore, if turbulent diffusion coefficients are converging, the ergodicity cannot be lost due to a limited sample size for a single trajectory. If diffusion coefficients do not converge (as it is often observed in numerical studies), then the diffusion process should be non-ergodic. Note that very often we can determine the ergodicity of the process from the analysis of behavior of ${\tau }_X^2$ (T, s) from s (see Section 7).

5. Lagrangian Trajectories

Two datasets of Lagrangian trajectories were used in this study to estimate the mean Lagrangian velocity and diffusion coefficients: surface drifters deployed in the Black Sea, and subsurface RAFOS floats launched in the California Current System.

The Black Sea drifters were Surface Velocity Program (SVP) drifters (NAVO, Washington, DC, USA; Marlin, Sevastopol, Ukraine; MetOcean, Dartmouth, NS, Canada) with a holey sock drogue at 15 m depth. The SVP drifters were satellite-tracked with 8–10 positions a day and geolocation accuracy of up to 1 km. For details on these data see [17,35]. The total dataset contains 44 trajectories collected in 1999–2003; we analyzed 15 trajectories (drifters #16330–#16337 and #33332–#33338) deployed in 2002–2003 because this time period contained most of the data (Figure 1a,b).

Drifter motions were highly coherent and determined by the basin and sub-basin scale structures of the Black Sea circulation (Figure 2): Rim current, gyres and semi-permanent anticyclonic eddies [36].

Isobaric RAFOS floats used in this study were mostly launched over the continental slope off Central California in 1992–2008. The durations of the trajectories vary from 10 days to ~3.3 years. RAFOS floats should be considered as single floats which cannot be analyzed by typical methods [18,31] because it is difficult to understand what should be taken as a statistical ensemble. Float positions have been determined with an absolute uncertainty not to exceed 10 km and about 1 km short-term relative error for individual floats [37]. After interpolating the trajectories with a 6-h step, diurnal tides and inertial motions have been filtered out.

In general, RAFOS floats moved poleward along the coast and thence westward; as the floats traveled, they were often trapped by different coherent structures, such as filaments, mesoscale and sub mesoscale eddies, etc. Here only two floats, #73 and #109 are used (see, Figure 3a,b). This shows [18] that float #73 was entrapped in two dynamic regimes: ballistic diffusion and sub-diffusion. Power exponents for these regimes are estimated below. The trajectory for float #109 suggested that it moved through similar regimes. Power exponents for float #109 are also estimated below.

6. Comparison with Other Methods

Calculations have demonstrated that the diffusion coefficients estimated by STM were at least one order smaller than estimates from the classical method from [10]. The difference between results obtained by STM and those from [10] and others using similar methods (see for example [35]) has a reasonable explanation: the difference is due to the way the observed trajectory was divided into its mean and diffusive parts. Except for the STM, the mean part of the trajectory is a Eulerian (or quasi-Eulerian) representation. For example, [10] used the quasi-Eulerian mean as the Lagrangian mean [1].

It was demonstrated [1] that for diffusion problems without any mean drift, approach [10] and wavelet decomposition give similar results. See, for example

Table 1. Comparison of diffusivity for Black Sea estimates using different mean circulation.

Method	K11(cm2/s) Zonal	K22(cm2/s) Meridional
Wavelet decomposition	7.5 × 106	2.6 × 106
Poulain et al. (2005) the mean is calculated by STM	4.5 × 106	1.4 × 106
Davis (1991) the mean is calculated by STM	6.9 × 106	2.4 × 106

. If we subtract a Eulerian mean (for example, taken from Demishev’s model [38] or the so-called pseudo-Eulerian mean [35]) from the total Lagrangian motions, results from [10] and the STM can be significantly different.

It is also not clear why in some cases the estimate of K_ij from [10] does not converge to a constant at quite long times (see Figure 6a–c from [1]). STM and the method from [10] are compared in greater detail in [1] (see Figures 5a–c and 14a–c in [1]).

7. Numerical Experiments with Black Sea Drifters and RAFOS Floats

7.1. Black Sea

The first group of numerical experiments was designed to test STM. It was applied to the seven and six trajectories shown in Figure 1a,b, respectively. It is expected that mean Lagrangian trajectories and turbulent diffusion coefficients calculated for individual Lagrangian trajectories should coincide for drifters deployed in the same area and at the same time.

Another problem here is how to estimate the anisotropy of the Black Sea turbulence and understand its spatial structure. The method developed by [15] is used to address this issue: (1) calculate the tensor of turbulent diffusion coefficients in cartesian coordinates; (2) find eigenvalues for this tensor; (3) calculate the ratio of the first and second eigenvalues. If this ratio is close to 1, the turbulent field is isotropic. Ratios larger or smaller than 1 indicate anisotropy of the turbulent field. This simple criterion should suffice for the analysis of drifter data.

The calculated elements of tensors K_ij and some mean characteristics for drifter #16330 are shown in Figure 4a,b, Figure 5a–d and Figure 6a–d; the results of calculations for all deployed drifters including coordinates of the middle points on drifter trajectories are given in

Table 2. Turbulent diffusion coefficients obtained by different methods for the first group of trajectories.

Drifter Number	Coefficients/106 (cm2/s) Calculated by the Wavelet Technique	Coefficients/106(cm2/s) Calculated from Equation (6)	Coordinates of the Middle Point (° E, ° N)
16330	K11$\approx 5$, K22$\approx 0.8,$K12$\approx$0.3	K11$\approx 0.5$, K22$\approx$0.5, K12$\approx$10−2	30.571; 43.0360
16331	K11$\approx 0.5$, K22$\approx$0.3, K12$\approx$0.1	K11$\approx$0.3, K22$\approx$0.1, K12$\approx$2 × 10−3	28.7060; 42.0040
16332	K11$\approx 6$,K22$\approx 2$, K12$\approx 0.2$	K11$\approx 0.2$, K22$\approx$0.1, K12$\approx$10−3	28.6940; 42.1040
16333	K11$\approx$2, K22$\approx 0.4$, K12$\approx$0.1	K11$\approx 1$, K22$\approx 0.8$, K12$\approx$6 × 10−3	29.2580; 42.2232
16334	K11$\approx 7$, K22$\approx 0.3$, K12$\approx$0.1	K11$\approx 2$, K22$\approx$0.4, K12$\approx$2 × 10−2	28.3948; 42.0032
16336	K11$\approx$2, K22$\approx$0.4, K12$\approx$0.1	K11$\approx 0.4$, K22$\approx$0.2, K12$\approx$10−2	29.994; 42.816

and

Table 3. Turbulent diffusion coefficients calculated for the second group of trajectories.

Drifter Number	Coefficients/106 (cm2/s) Calculated by the Wavelet Technique	Coefficients/106(cm2/s) Calculated from Equation (6)	Coordinates of the Middle Point (° E, ° N)
33347	${\mathrm{K}}_{11} \approx$2.0, K22$\approx 0.8$,${\mathrm{K}}_{12} \approx$0.1	${\mathrm{K}}_{11} \approx$0.6, K22$\approx$0.2, ${\mathrm{K}}_{12} \approx$10−3	34.233; 43.981
33348	-	${\mathrm{K}}_{11} \approx$10, K22$\approx$2, ${\mathrm{K}}_{12} \approx$0.3	37.166; 44.794
33349	${\mathrm{K}}_{11} \approx$2.0, K22$\approx 0.3$,${\mathrm{K}}_{12} \approx$0.3	${\mathrm{K}}_{11} \approx$1.0, K22$\approx$0.6, ${\mathrm{K}}_{12} \approx$10−3	34.4132; 44.301
33350	${\mathrm{K}}_{11} \approx$0.8, K22$\approx$0.5, ${\mathrm{K}}_{12} \approx$0.4	${\mathrm{K}}_{11} \approx$0.4, K22$\approx$0.2, ${\mathrm{K}}_{12} \approx$2 × 10−3	34,4132; 44,1077
33351	${\mathrm{K}}_{11} \approx$1.0, K22$\approx$0.5, ${\mathrm{K}}_{12} \approx$0.5	${\mathrm{K}}_{11} \approx$0.4, K22$\approx$0.1, ${\mathrm{K}}_{12} \approx$0.5 × 10−3	35.828; 44.47
33352	${\mathrm{K}}_{11} \approx$1.0, K22$\approx$0.8, ${\mathrm{K}}_{12} \approx$4 × 10−2	${\mathrm{K}}_{11} \approx$0.4, K22$\approx$0.2, ${\mathrm{K}}_{12} \approx$10−3	33.95; 44.38

.

Figure 4a,b compares the mean Lagrangian and observed trajectories of drifters #16330 and #33347. It is obvious from this figure that the residual (a difference between the observed and mean Lagrangian trajectories) for these trajectories is quite small. This demonstrates a large degree of smoothing of Lagrangian data.

The very short trajectory of drifter #16335 was excluded from calculations. We see from Table 2 that the drifters can be combined into two groups: (#16330, #16333, #16336) and (#16331, #16332, #16334). This gives the coefficients of turbulent diffusion in two areas:

$$K = \left(\begin{array}{cc}4.5& 0.15\\ 0.15& 0.5\end{array}\right) \times {10}^6 \mathrm{and} {\lambda }_1/{\lambda }_2 = 0.5/4.5\approx 0.11 \mathrm{for} \mathrm{the} \mathrm{first} \mathrm{group} \mathrm{of} \mathrm{drifters} \mathrm{in} \mathrm{the} \mathrm{area} \mathrm{near} 29.148° \mathrm{E}, 41.763° \mathrm{N},$$

(25)

$$K = \left(\begin{array}{cc}2.9& 1.3\\ 1.3& 0.8\end{array}\right) \times {10}^6 \mathrm{and} {\lambda }_1/{\lambda }_2 = 0.18/3.52\approx 0.05 \mathrm{for} \mathrm{the} \mathrm{second} \mathrm{group} \mathrm{of} \mathrm{drifters} \mathrm{in} \mathrm{area} \mathrm{near} 30.243° \mathrm{E}; 43.445° \mathrm{N}.$$

(26)

Estimates (25) and (26) indicate that the turbulence in the western part of the Black Sea was strongly anisotropic. However, it is stronger off the coast of Turkey ($\lambda$₁/$\lambda$₂$\approx 0.05$) than near Crimea ($\lambda$₁/$\lambda$₂$\approx 0.11$). The analysis of the western Black Sea circulation show that the first eigenvector is directed perpendicular to the direction of the Rim current.

A second group of numerical experiments was applied to estimate the turbulent diffusion tensor near Crimea (in the area of the Kerch strait). Six drifters (#33347–#33352) were deployed in this area and moved westward. Trajectories of drifters #33348 and #33351 were very short and that resulted in reduced accuracy of estimates for high order moments of ${\tau }_X^2$ (T, s). Drifter #33348 moved only 12 days and was excluded from the analysis.

Estimated eddy turbulent diffusion coefficients are plotted in Figure 6a–d and summarized in Table 3. Calculations were similar to those done for Lagrangian drifters in Table 2. Table 3 shows that drifters #33347, #33349, #33350, and #33352 can be grouped. This gives the coefficients of turbulent diffusion in two different areas as

$$K = \left(\begin{array}{cc}6.8& 0.32\\ 0.32& 1.0\end{array}\right) \times {10}^6 \mathrm{and} {\lambda }_1/{\lambda }_2 = 0.98/6.8\approx 0.14 \mathrm{for} \mathrm{the} \mathrm{first} \mathrm{group} \mathrm{of} \mathrm{drifters} \mathrm{in} \mathrm{the} \mathrm{area} \mathrm{near} 34.546° \mathrm{E} 44.150° \mathrm{N},$$

(27)

$$K = \left(\begin{array}{cc}1.0& 0.5\\ 0.5& 0.5\end{array}\right) \times {10}^6 \mathrm{and} {\lambda }_1/{\lambda }_2 = 0.19/1.3\approx 0.15 \mathrm{for} \mathrm{the} \sec \mathrm{ond} \mathrm{group} \mathrm{of} \mathrm{drifters} \mathrm{in} \mathrm{the} \mathrm{area} \mathrm{near} 37.111° \mathrm{E} 44.592° \mathrm{N}.$$

(28)

Based on the $\lambda$₁/$\lambda$₂ ratios, the turbulence generated in these areas by large eddies is also anisotropic. Note that trajectories of all drifters deployed in the Black Sea were quasi-ergodic as demonstrated by Figure 5b and Figure 6b and follows from [20].

There is no numerical bias due to the limited number of drifters/floats because we work only with individual trajectories. However, the estimates of the turbulent diffusion coefficients may differ for different sample lengths of a drifter/float trajectory. As seen from an example shown in Figure 5c,d, turbulent diffusion coefficients tend to be constants for t > 100 days. Estimates of coefficients calculated for t = 40 and for t = 150 days may differ by as much as five times. From our point of view, the length of a trajectory longer than 100–150 days is quite satisfactory to obtain estimates which are non-sensitive to the length of drifter/RAFOS float trajectories [39].

7.2. California Current System

The trajectories of RAFOS floats #73 and #109 which were deployed in the California Current system are examined next (Figure 7a–d). Results obtained by [18] have shown that the first float demonstrated non-Gaussian anisotropic sub-diffusion behavior for t > ${\tau }_0$ and non-Gaussian isotropic ballistic diffusion for t < ${\tau }_0$ (see Figures 16d and 17d from [18,40]), where ${\tau }_0$ ~20–30 days. Visual analysis demonstrates that the floats were trapped by eddy-like structures.

Calculations show [18] that for anisotropic sub-diffusion, elements of the dispersion tensor are $\sigma$₁₁~t^q and $\sigma$₂₂~t^p and q = 0.78 (0.38) and p = 0.15 (0.05) for float #73 (#109). The meridional diffusion can be neglected and the equations for particle dispersion rewritten as $\sigma$_eff~t^q and then the estimate of the power exponent of q [41]:

$$\mathrm{q} = \mathrm{t}/\langle \delta X^2\rangle \times \mathrm{d} \langle \delta X^2\rangle /\mathrm{dt},$$

(29)

As expected for sub-diffusion, this method gives constant power exponents. However, for the transition between ballistic (or normal) and sub-diffusion regimes, exponents are not constant.

To introduce additional averaging, each float trajectory was divided into 15 sub- trajectories, each approximately 25.5 days, and this ensemble was used for additional averaging <<…>>. The results were not too sensitive to an increase in the length of each sub trajectory. Therefore, no additional mathematical procedures were used to increase the stability of the results.

RAFOS float #73. The results (Figure 8a) showed two different dynamical regimes: normal diffusion (not ballistic as in [18]) for t < 10 days ($\alpha \approx 1)$ and sub-diffusion for t > 20 days ($\alpha \approx 0.78).$ Our estimates also demonstrated that motions of float #73 were neither ergodic (EB$\approx$0.25) nor Gaussian (H$\approx$10²) in both regimes.

RAFOS float #109. This float was not included in [18] because it did not complete its mission until 2006. Since float #109 trajectory was longer than #73, the length of sub-trajectories was increased to 50 days. The normal diffusion regime was observed for approximately 30 days but starting on day 45 it was replaced by sub-diffusion with $\alpha \approx$ 0.63 (see Figure 8b). The float #109 trajectory was neither ergodic (EB $\approx$ 0.3) nor Gaussian (H $\approx$ 10²).

Therefore, two RAFOS floats demonstrated non-Gaussian and non-ergodic behavior after being trapped by mesoscale eddy-like structures. Both the floats indicated normal diffusion followed by sub-diffusion. Similar power exponents were estimated for both cases.

8. Conclusions

The STM was applied to three groups of Lagrangian data collected from single drifters/floats deployed in the Black Sea and California Current system. The STM technique replaces averaging over a statistical ensemble by averaging along a single Lagrangian trajectory for ergodic (normal diffusion) and weakly non-ergodic (sub-diffusion) flows. Criteria to estimate the degree of non-ergodicity and non-Gaussivity for flows are formulated. The technique provides estimates of power exponents for the normal diffusion and sub-diffusion processes.

STM results demonstrated that estimates of reconstructed diffusion coefficients strongly depended on how the trajectory was divided into its mean Lagrangian and diffusion parts. An incorrect choice of the mean Lagrangian trajectory may result in erroneous estimates of turbulent diffusion coefficients.

Black Sea turbulent diffusion was strongly anisotropic, and the strength of anisotropy varied along the direction of the Rim current. On average, the turbulent diffusion coefficients generated by intense turbulent eddies and waves were approximately 10⁵–10⁶ cm²/s. Estimates obtained by [35] for the same Lagrangian data were considerably higher, about 10⁷–10⁸ cm²/s. This difference is attributed to [35] use of the Davis (1991) technique to divide the flow into quasi-Eulerian mean and diffusion parts.

Normal diffusion and sub-diffusion were studied using two isobaric subsurface RAFOS floats in the California Current system. Here a method for estimating power exponents was developed which was constant for the normal diffusion and sub-diffusion regimes but varied in the transition between these two regimes.

The STM method assumed that Lagrangian particle movement can be divided into two parts: mean Lagrangian drift and turbulent diffusion. More accurately, these diffusive motions can be described as a process with memory [26]. For example, the Mori-Zwanzig formalism seems to be very useful in such situations. From the theoretical point of view, this formalism yields equations which correctly describe particle dynamics and asymptotes for particle behavior when t$\to \infty .$

The STM approach for the analysis of single drifter/float trajectories is computationally simple and has the following advantages:

The approach developed in [1] and the present paper seem to be relatively universal and generalizable for a multi-dimensional diffusivity, as, for example, in multidimensional atmospheric models [42,43]. Applying the wavelet decomposition allows for dividing single Lagrangian trajectories into their mean and diffusive parts. The former represents motions with degrees of freedom larger than one or two.