2.1. Algorithm

Let us consider the LM of the field $\phi (\boldsymbol {x},t)$ from $t_0$ to $t_0+T$. Before explaining the steps of GBLA, we clarify some notations and definitions. We denote the Eulerian mean of $\phi (\boldsymbol {x},t)$ over the period of $t_0$ to $t$ with

(2.1)\begin{equation} \bar{\phi}(\boldsymbol{x},t_0,t) = \frac{1}{t-t_0} \int_{t_0}^{t} \phi(\boldsymbol{x},\tau) \,{\rm d} \tau. \end{equation}

The LM of $\phi$ is the time average of this quantity along the trajectory, that is, at fixed particle label instead of fixed spatial position. There are different ways of labelling particles; for instance, the particle position at the beginning or the end of the averaging interval may be selected as labels. Denoting this label by $\boldsymbol {x}$, we can map it to the particle position at time $t$ (denoted by $\boldsymbol {X}$) using the particle displacement field $\boldsymbol {{\xi }}(\boldsymbol {x},t)$,

(2.2)\begin{equation} \boldsymbol{X} = \boldsymbol{{\Xi}}(\boldsymbol{x},t) = \boldsymbol{x}+\boldsymbol{{\xi}}(\boldsymbol{x},t). \end{equation}

Using (2.1) and (2.2), the LM operator is mathematically described as

(2.3)\begin{equation} \overline{\phi}^{{\hspace{0.5mm}\small L}} (\boldsymbol{x},t_0,t) = \frac{1}{t-t_0} \int_{t_0}^{t} \phi(\boldsymbol{x}+\boldsymbol{{\xi}}(\boldsymbol{x},\tau),\tau) \,{\rm d} \tau. \end{equation}

Different choices of $\boldsymbol {x}$ lead to different displacement fields $\boldsymbol {{\xi }}$. Andrews & McIntyre (Reference Andrews and McIntyre1978) defined the ‘generalised Lagrangian mean (GLM)’ by requiring

(2.4)\begin{equation} \bar{\boldsymbol{{\xi}}}(\boldsymbol{x},t) = 0. \end{equation}

In other words, $\boldsymbol {x}$ is now the mean position of the particle from $t_0$ to $t$. If the displacement field is small (compared with the changes in mean position), (2.4) gives the LM some Eulerian-like character as $\boldsymbol {x}$ is more than just a label and represents the spatial location of the trajectory for which the mean is calculated. As it will be shown in the examples of § 3.1, this choice of label is more effective at removing the fast time scale in the LM estimate provided that the mapping (2.2) does not have a fast time dependence. If other labels such as instantaneous position at the end or middle of the averaging interval are chosen, the filtered signal may still contain some fast oscillations or get distorted (see the numerical example of § 3.1). Moreover, the GLM definition holds properties that lead to significant analytical simplifications in modelling two time scale geophysical flows (Andrews & McIntyre Reference Andrews and McIntyre1978; Bühler & McIntyre Reference Bühler and McIntyre1998; Xie & Vanneste Reference Xie and Vanneste2015; Gilbert & Vanneste Reference Gilbert and Vanneste2018). Throughout this paper, we use (2.3) in conjunction with (2.4) to express the LM in terms of the mean position. For $t<(t_0+T)$, (2.1) and (2.3) represent partial means, and for $t = t_0+T$ the complete mean for the intended averaging interval. Considering that in this section we focus only on one averaging interval with a fixed starting point, we drop $t_0$ in the arguments of the mean operators to lighten the notation.

The core idea behind GBLA is to calculate partial LMs recursively in time. For instance, we assume that the partial LM from $t_0$ to $t_k$ ($t_0< t_k< t_0+T$) for the particles that reach the grid points at time $t_k$ is already computed. Through the steps laid out below, the LM is extended to $t_{k+1}$ and calculated for a new set of particles that reach the grid points at time $t_{k+1}$ (using the partial averages of the previous step). In this approach, at each time step a different set of trajectories is considered (the set that ends up at grid points at that time step), which is in contrast with particle tracking methods where one trajectory is followed all along over the entire averaging interval. The incremental extension of partial LMs continues to the end of the averaging interval $t_0+T$. The completed LM field is then assigned to time $t_0+T/2$. Note that in our demonstration we focus only on one averaging interval, but this process can be repeated for new intervals, and averaging procedure can be carried out simultaneously. For example, the next interval may start from $t_m$ and end at $t_{m+N}$ and so on. If the beginning of a later interval exceeds $t_0+T$, the memory allocated to the first interval (from $t_0$ to $t_0+T$) can be freed in a systematic way and used for the new intervals.

Equipped with the above definitions, we present the steps to compute $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {X}_i,t_{k+1})$, $t_{k+1})$ from $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\{\boldsymbol {X}_{\alpha }\},t_k), t_k)$, where $\boldsymbol {X}_i$ is an arbitrary grid point in space and $\{\boldsymbol {X}_{\alpha }\}$ is a set consisting of $\boldsymbol {X}_i$ and its neighbours. In other words, our steps show how to calculate the partial LM (from $t_0$ to $t_{k+1}$) for the particle that passes through $\boldsymbol {X}_i$ at $t_{k+1}$ by knowing the LM (from $t_0$ to $t_k$) for the particles that pass through $\boldsymbol {X}_i$ and its neighbouring grid points at $t_{k}$. To simplify the notation we use one spatial index $i$, which can be viewed as a multielement vector or replaced by several indices in higher than one dimension. Considering that $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {X}_i,t_k),t_k)$ is the partial LM for the particle that reaches $X_i$ at $t_{k}$ and $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {X}_i,t_{k+1}),t_{k+1})$ is the partial LM for the particle that reaches $X_i$ at $t_{k+1}$, these two LMs belong to two separate trajectories (marked in red and blue in figure 1).

1. Step 1 Calculate $\phi (\boldsymbol {X}_i, t_{k+1})$ and the velocity $\boldsymbol {u}(\boldsymbol {X}_i, t_{k+1})$ by integrating the governing equations from $t_k$ to $t_{k+1}$. Note that a finer time step may be used for the time integration of the governing equations than that used for averaging; there might be several simulation time steps between $t_k$ and $t_{k+1}$.

2. Step 2 Consider the particle that goes through $\boldsymbol {X}_i$ at time $t_{k+1}$ (red trajectory in figure 1). Approximate the position of this particle at time $t_k$ by advecting it backward in time,

   (2.5)\begin{equation} \boldsymbol{Y}_i = \boldsymbol{X}_i - (t_{k+1}-t_k) \boldsymbol{u}(\boldsymbol{X}_i, t_{k+1}). \end{equation}
   This point is marked in green in figure 1 and splits the red trajectory into two parts.

3. Step 3 Denote the partial LM of $\phi$ between $t_0$ and $t_k$ for the red trajectory by $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {Y}_i,t_k),t_k)$. Then, find its value by interpolation using the LM for the trajectories that go through the neighbouring grid points (marked by blue curves). For example, $\boldsymbol {Y}_i$ lies between $\{\boldsymbol {X}_{\alpha }\}$ for ${\alpha } = i, j, l$ and $p$ in figure 1. Therefore, the four values of $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\{\boldsymbol {X}_{\alpha }\},t_k),t_k)$ may be used to interpolate $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {Y}_i,t_k),t_k)$. Note that $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {Y}_i,t_k),t_k)$ and $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {X}_i,t_k),t_k)$ are both averages over $[t_0,\ t_k]$ but for different trajectories (highlighted in orange and light blue, respectively).

4. Step 4 Selecting the mean position of the red trajectory from $t_0$ to $t_{k+1}$ as a label for the particle that passes through $\boldsymbol {X}_i$ at $t_{k+1}$ and denoting it by $\boldsymbol {x}_i$, use the definition (2.3) to compute

   (2.6)\begin{align} &\overline{\phi}^{{\hspace{0.5mm}\small L}} (\boldsymbol{{\Xi}}^{{-}1}(\boldsymbol{X}_i,t_{k+1}),t_{k+1}) = \frac{1}{t_{k+1}-t_0} \int_{t_0}^{t_{k+1}} \phi(\boldsymbol{x}_i + \boldsymbol{{\xi}}(\boldsymbol{x}_i,\eta),\eta) \,{\rm d} \eta\nonumber\\ & = \frac{t_k-t_0}{t_{k+1} - t_0} \ \frac{1}{t_{k} - t_0} \int_{t_0}^{t_k} \phi(\boldsymbol{x}_i + \boldsymbol{{\xi}}(\boldsymbol{x}_i,\eta),\eta) \,{\rm d} \eta+ \frac{1}{t_{k+1} - t_0} \int_{t_k}^{t_{k+1}} \phi(\boldsymbol{x}_i + \boldsymbol{{\xi}}(\boldsymbol{x}_i,\eta),\eta) \,{\rm d} \eta\nonumber\\ & = \frac{t_k-t_0}{t_{k+1} - t_0}\ \overline{\phi}^{{\hspace{0.5mm}\small L}} (\boldsymbol{{\Xi}}^{{-}1}(\boldsymbol{Y}_i,t_k),t_k) +\frac{1}{t_{k+1} - t_0} \int_{t_k}^{t_{k+1}} \phi(\boldsymbol{x}_i + \boldsymbol{{\xi}}(\boldsymbol{x}_i,\eta),\eta) \,{\rm d} \eta. \end{align}
   The first term on the right-hand side is already computed in the previous step. The second term, which is the weighted mean of the red trajectory from $t_k$ to $t_{k+1}$, can be approximated by
   (2.7)\begin{equation} \int_{t_k}^{t_{k+1}} \phi(\boldsymbol{x}_i + \boldsymbol{{\xi}}(\boldsymbol{x}_i,\eta),\eta) \,{\rm d} \eta \approx \phi(\boldsymbol{X}_i,t_{k+1})\ (t_{k+1}-t_k). \end{equation}
   A better approximation may be derived by the trapezoidal rule
   (2.8)\begin{equation} \int_{t_k}^{t_{k+1}} \phi(\boldsymbol{x}_i + \boldsymbol{{\xi}}(\boldsymbol{x}_i,\eta),\eta) \,{\rm d} \eta \approx \tfrac{1}{2}\left[ \phi(\boldsymbol{Y}_i,t_{k}) +\phi(\boldsymbol{X}_i,t_{k+1}) \right] (t_{k+1}-t_k). \end{equation}
   To calculate $\phi (\boldsymbol {Y}_i,t_{k})$, another interpolation is required. Hence, depending on the desired accuracy and computational cost, (2.8) may or may not be preferred to (2.7).
   Figure 1.

   Schematic view of particle trajectories that pass through the grid points $\boldsymbol {X}_i$, $\boldsymbol {X}_j$, $\boldsymbol {X}_p$ and $\boldsymbol {X}_l$ that are marked by black. The trajectories that reach these grid points at $t_k$ are marked in blue and the trajectory that reaches the grid point at $t_{k+1}$ in red. For $\boldsymbol {X}_i$, the times associated with the particle position of red and blue trajectories are shown as well.

   

   So far the calculation of the LM field in (2.6) is complete, but it is given for the end position of each particle $\boldsymbol {X}_i$. If the definition of GLM with the condition (2.4) is desired, the LM should be assigned to the mean position of particles. In other words, we need to calculate $\boldsymbol {x}_i = \boldsymbol {{\Xi }}^{-1}(\boldsymbol {X}_i,t_{k+1})$, which maps the particle's end position $\boldsymbol {X}_i$ to its mean from $t_0$ to $t_{k+1}$. To calculate the mean position, we update it in a similar manner to (2.6) at each time step via step 5.

5. Step 5 Assuming $\boldsymbol {{\Xi }}^{-1}(\boldsymbol {X},t_{k})$ is known for all grid points from the previous time step, compute the mean position of the red trajectory from $t_0$ to $t_{k+1}$,

   (2.9)\begin{align} \boldsymbol{{\Xi}}^{{-}1}(\boldsymbol{X}_i,t_{k+1}) &= \frac{1}{t_{k+1}-t_0} \int_{t_0}^{t_{k+1}} \boldsymbol{x}_i + \boldsymbol{{\xi}}(\boldsymbol{x}_i,\eta) \,{\rm d} \eta\nonumber\\ & = \frac{t_k-t_0}{t_{k+1} - t_0}\ \boldsymbol{{\Xi}}^{{-}1}(\boldsymbol{Y}_i,t_{k}) +\frac{1}{t_{k+1} - t_0} \frac{1}{2}\left( \boldsymbol{Y}_i + \boldsymbol{X}_i \right) (t_{k+1}-t_k), \end{align}
   where $\boldsymbol {{\Xi }}^{-1}(\boldsymbol {Y}_i,t_{k})$ is interpolated using the neighbouring grid points. Note that (2.9) is simply a particular case of (2.6) in conjunction with (2.8) for $\phi (\boldsymbol {x},t) = \boldsymbol {x}$.

The above steps are carried out for all the grid points and iteratively in time until $t_{k+1} = t_0+T$. The completed LM is now given on scattered mean positions $\boldsymbol {{\Xi }}^{-1}(\boldsymbol {X}_i,t_{N})$ and can be interpolated back on the computational grid points.

2.2. Parameter set-up, convergence and numerical error

The GBLA method may be implemented on temporally and spatially decimated fields. To distinguish between the two, we denote the (iterative) averaging time step with $\Delta t = t_{k+1} - t_k$ and the simulation time step (used to solve the governing equations) with $\delta t$ throughout this paper. Similarly, $\Delta \boldsymbol {x}$ denotes the spatial resolution of GBLA, which could be coarser than the one used in simulations (i.e. weather or climate models). It is important to keep $\Delta t$ and $\Delta \boldsymbol {x}$ smaller than the characteristic time and length scales of the fast phenomenon that is going to be filtered. To guarantee the stability of backward advection in (2.5) the Courant–Friedrichs–Lewy condition $\boldsymbol {u} \Delta t < \Delta \boldsymbol {x}$ should be satisfied.

Just like any numerical algorithm, GBLA is susceptible to numerical errors. Denoting the local error at each iteration (from $t_k$ to $t_{k+1}$) by $e_k$, there are three main contributors to $e_k$.

1. (i) Finite differentiation used to advect particles backward in (2.5). Considering the Taylor expansion of particle position on a trajectory, the difference between the exact and approximated particle position at time $t_k$ is $O(\Delta t^2)$.

2. (ii) Interpolation of LM at time $t_k$. In the step 3 of the algorithm, we interpolate $\overline {\phi }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {Y}_i,t_k),t_k)$, leading to an error that depends on the selected interpolation scheme. Considering polynomial interpolation of degree $n$ on a grid with spacing of $\Delta x$, one can show that the interpolation error is $O((\Delta x)^{n+1})$. For example, the linear interpolation error is $O((\Delta x)^2)$ and the cubic interpolation error $O((\Delta x)^3)$. Note that the order of error is independent of spatial dimension.

3. (iii) Approximation of the last integral in (2.6) either by (2.7) or (2.8). The error of the former is $O(\Delta t)$ and the latter $O(\Delta t^3 (\Delta x)^{n+1})$ considering the interpolation required to calculate $\phi (\boldsymbol {Y}_i,t_{k})$.

Considering that there are linear leading terms in polynomial interpolation, the error in (i) directly translates into the total local error $e_k$. As a result, we can write

(2.10)\begin{equation} e_k = O(\Delta t^2) + O((\Delta x)^{n+1}) + O(\Delta t^s), \end{equation}

where $s =1$ or $3$ if (2.7) or (2.8) are used, respectively. Hence, depending on the value of $s$, the order of local error in $\Delta t$ may be determined by (i) or (ii). The calculation of global error, which is the accumulation of local errors, is more complicated. It cannot be derived simply by multiplying $e_k$ with the number of time steps $N = T/\Delta t$, because in the interpolation of the step 3 the neighbouring points contain the accumulated interpolation errors from the previous iterations. As a result, different sources of error become intertwined. Comparatively, the particle tracking approach does not have this type of error accumulation, because it uses the instantaneous values at neighbouring points for each time step. Hence, GBLA might be more affected by interpolation errors than the particle tracking methods. We leave the analytical expression of global error for a future study, but in § 3.1 and § 3.2.1 we calculate it numerically to investigate the convergence of GBLA.

3.1. Illustrative one-dimensional flows

We consider the following example with a velocity that consists of a fast and slow component:

(3.1a–c)\begin{equation} u = u_{f} + u_{s}, \quad u_{f} = \epsilon \sin(20x - 20t), \quad u_{s} = \mathrm{e}^{{-}a(x-{\rm \pi})^2}, \end{equation}

where $a = 15$ for $x<{\rm \pi}$ and $a = 1$ for $x\geq {\rm \pi}$. Figure 2(a) displays the LM velocity of this flow regarded as a function of mean position using GBLA with cubic interpolation (shown in magenta) and particle tracking (shown in black). For both methods, $\Delta x = 4.2\times 10^{-3}$, $\Delta t = 2.6\times 10^{-3}$ and the averaging period $T = 2{\rm \pi} /3$, which is longer than the period of the fast signal $2{\rm \pi} /20$. The GBLA and particle tracking results (magenta and black curves) perfectly lie on top of each other, which confirms the validity of GBLA. To better quantify the accuracy and convergence of GBLA, in figure 2(b) we coarsen $\Delta x$ and calculate the following error norm (representing the global error for the averaging period of $T$):

(3.2)\begin{equation} \Vert e \Vert=\sqrt{\frac{1}{N} \sum (\overline{u}^{{\hspace{0.5mm}\small L}} _{{GBLA}} - \overline{u}^{{\hspace{0.5mm}\small L}} _{{exact}})^2}, \end{equation}

where the summation is carried out over all spatial grid points, and $N$ is the number of spatial grid points. Here $\overline {u}^{{\hspace {0.5mm}\small L}} _{{GBLA}}$ is the LM calculated using GBLA and $\overline {u}^{{\hspace {0.5mm}\small L}} _{{exact}}$ is the exact value of LM which is approximated by particle tracking at much finer spatial and temporal resolutions of $\Delta x = 2 \times 10^{-3}$ and $\Delta t = 2.6 \times 10^{-3}$. Two different interpolation schemes are used in computing GBLA for figure 2(b): linear (black dots) and cubic (red dots). The overall message of this figure is that GBLA with different interpolation schemes converges to the exact solution if $\Delta x$ is lowered. At smaller $\Delta x$, GBLA's error with cubic interpolation has a slope close to 2, and with linear interpolation a slope slightly shallower than 1. At larger $\Delta x$ other sources of error contribute to $\Vert e \Vert$ leading to shallower slopes.

Figure 2.

(a) The LM velocity of (3.1a–c) calculated for $\epsilon = 0.06$ from $t=0$ to $t=2{\rm \pi} /3$ for $\Delta x = 4.2\times 10^{-3}$ and $\Delta t = 2.6\times 10^{-3}$ using: GBLA (with cubic interpolation) as a function of mean position (magenta); particle tracking as a function of mean position (black); particle tracking as a function of position at the middle of the averaging interval (blue); and particle tracking as a function of position at the end of the averaging interval (green). (b) The error norm $\Vert e \Vert$ of the GBLA scheme using linear interpolation (black dots) and cubic interpolation (red dots). The slope of blue guide line is 1 and that of red line is 2.

As mentioned in the previous section, we choose the definition of GLM theory by requiring (2.4) to hold and using mean position as a label for LM. We use the above example to clarify the importance of this choice. In figure 2(a), in addition to expressing the LM as a function of mean position, we show the LM as a function of particle position at the middle of the averaging interval (blue curve) and as a function of particle position at the end of the averaging interval (green curve). We have used particle tracking methods to calculate these results. The black and blue curves preserve the overall shape of the slow signal, but the green curve fails to do so, which shows that the LM in terms of end position may distort the mean-flow signal. The black and blue curves have differences as well, which are more visible at the steep velocity gradient before $x={\rm \pi}$. To understand the reason behind this discrepancy, consider a particle that travels through the relatively flat part of velocity profile in the first half of the averaging interval and through the steep velocity jump in the second half. As a result, the average speed of this particle is small in the first half and large in the second half, leading to a relatively small change of position in the first half compared with the second half. Therefore, the mean position of this particle differs from its position at the middle of the averaging interval (which is closer to the position at the beginning of the interval).

In figure 3 we increased the wave amplitude to $\epsilon = 0.15$ to illustrate how the GLM definition of LM can be more effective in removing the fast time scale (all other parameters are the same as those used for figure 2a). The LM at different times is shown in terms of the mean position calculated by GBLA (magenta curve) and in terms of the position at the middle of the averaging interval calculated by particle tracking (blue curve). There is still a noticeable trace of fast signal in the blue curves, which propagates rightward in time. This is because the particle position at the middle of the averaging period oscillates with the fast time scale, whereas the mean position is less affected by the these oscillations. This example is deliberately chosen to amplify the difference for various choices of spatial labels. For smaller-amplitude waves and milder background velocity gradients the magenta and blue curves converge to each other. Applying Lagrangian filtering to oceanic flows, many authors have used the final position as a label for LM (Shakespeare & Hogg Reference Shakespeare and Hogg2017, Reference Shakespeare and Hogg2019) or the position at the middle of the averaging interval (Shakespeare et al. Reference Shakespeare, Gibson, Hogg, Keating, Bachman and Velzeboer2021). If a more accurate removal of fast time scale is desired, expressing LM in terms of mean position should be considered instead.

Figure 3.

The LM velocity of (3.1a–c) over the period of $2{\rm \pi} /3$, calculated for $\epsilon = 0.15$ at different times using: GBLA in terms of mean position (magenta) and particle tracking in terms of position at the middle of the averaging interval (blue).

3.2. Using a materially conserved quantity to investigate the validity of GBLA

In this section, we consider the materially conserved quantities – which are constant on each trajectory – to study the accuracy of GBLA. For these quantities, the LM expressed in terms of final position should be equal to their instantaneous values at the end of the averaging interval. We employ this fact to see how reliable GBLA is. We remind the reader that this is merely for studying the validity of GBLA; as discussed before using the GLM definition and expressing the LM in terms of mean position instead of final position is preferred when it comes to the applications of Lagrangian averaging. As a result, step 5 of § 2 is skipped in the calculations of this part.

3.2.1. Conservation of vorticity in two-dimensional incompressible flow

For two-dimensional incompressible inviscid flows, vorticity is conserved on particle trajectories,

(3.3a,b)\begin{equation} \frac{\partial \zeta}{\partial t} + u \frac{\partial \zeta}{\partial x} + v \frac{\partial v}{\partial y} = 0, \quad \zeta = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}. \end{equation}

As a result, $\zeta (\boldsymbol {X},t) = \overline {\zeta }^{{\hspace {0.5mm}\small L}} (\boldsymbol {{\Xi }}^{-1}(\boldsymbol {X},t),t)$, where $\boldsymbol {X} = (x,y)$. Therefore, to examine GBLA we derive the LM vorticity (expressed in terms of the final particle position) and compare it with the instantaneous vorticity field at the end of the averaging interval. We initialise the numerical simulation with two like-signed vortices,

(3.4) $$\begin{gather} \zeta(x,y,t=0) = \exp({-(x-{\rm \pi}+0.1)^2 -(y-{\rm \pi}+{\rm \pi}/3)^2 })\nonumber\\ + \exp({-(x-{\rm \pi}-0.1)^2 - (y-{\rm \pi}-{\rm \pi}/3)^2}). \end{gather}$$

We employ a standard pseudo-spectral code to solve (3.3a,b) on a $2{\rm \pi} \times 2{\rm \pi}$ domain with periodic boundary conditions using $256 \times 256$ grids points. The time step for implementing the averaging process is 10 times coarser than the simulation time step, $\delta t = 0.002$. Figure 4 portrays the LM vorticity calculated by different methods next to the instantaneous vorticity at the end of the averaging interval for two different averaging periods. For the shorter averaging interval, all methods lead to very similar results. However, for the longer interval, GBLA with linear interpolation is less accurate than the others. The comparison between panels (c) and (h), and (d) and (i), shows that the use (2.8) instead of (2.7) has not made any discernible difference, because the interpolation error exceeds other types of error due to the relatively small time step.

Figure 4.

Instantaneous and LM vorticity field for two like-signed vortices. Panels (a–e) correspond to fields at $t=10$ and ( f–j) to $t=30$. (a, f) Instantaneous vorticity; (b,g) GBLA using linear interpolation and (2.7); (c,h) GBLA using cubic interpolation and (2.7); (d,i) GBLA using cubic interpolation and (2.8); (c,h) LM by particle tracking and using cubic interpolation.

To have a more quantitative comparison between different methods, we calculate the following error norm (a two-dimensional version of the norm used in (3.2)):

(3.5)\begin{equation} \Vert e \Vert=\sqrt{\frac{1}{N_x N_y} \sum_i (\overline{\zeta}^{{\hspace{0.5mm}\small L}} (\boldsymbol{{\Xi}}^{{-}1}(\boldsymbol{X}_i,T),T)-\zeta(\boldsymbol{X}_i,T))^2 }, \end{equation}

where the summation is carried out over all spatial grid points, and $N_x$ and $N_y$ are the number of grid points in each dimension. Here $\zeta (\boldsymbol {X}_i,T)$ is an approximation for the exact LM, which is obtained by numerically solving (3.3a,b) at a high resolution of $384 \times 384$ and small time step of $\delta t = 0.002$. Here $\overline {\zeta }^{{\hspace {0.5mm}\small L}}$ is calculated using different Lagrangian averaging methods with temporal and spatial distancing that are lowered to investigate their convergence and accuracy.

In figure 5(a) we set the averaging time step to $\Delta t = 0.02 = 10\, \delta t$, which is small enough to highlight the interpolation error (see § 2.2 for discussion on different sources of error). The GBLA method with linear interpolation and using (2.7) (red dots) underperforms GBLA and particle tracking with cubic interpolation (black squares and blue circles, respectively). However, after using (2.8) even GBLA with linear interpolation (magenta triangles) yields accurate results within the same ‘ballpark’ as cubic interpolation and particle tracking approaches. With cubic interpolation, the convergence of GBLA is $O((\Delta x)^2)$, and that of particle tracking is $O(\Delta x)$. This is similar to the spatial convergence of error observed in figure 2(b) for the one-dimensional synthetic flow. For large $\Delta x$, GBLA with cubic interpolation has higher errors than its particle tracking counterpart, but by lowering $\Delta x$ both methods converge to the same result.

Figure 5.

(a) Global error using different spatial resolution. (b) Global error using different averaging time steps. Error norm is calculated for different schemes namely: GBLA with linear interpolation and using (2.7) (red dots); GBLA with cubic interpolation and using (2.7) (black squares); particle tracking with cubic interpolation (blue circles); and GBLA with linear interpolation and using (2.8) (magenta triangles). The blue dashed lines have the slope of 1, and the black line has the slope of 2.

In figure 5(b) the convergence of different methods with respect to $\Delta t$ is investigated. For all the calculations of this plot $\Delta x = \Delta y = 2{\rm \pi} /384$. Similar to the figure 5(a), the error of GBLA with linear interpolation and using (2.7) (red dots) is almost two orders of magnitude higher than the other methods. Use of (2.8) reduces the error of linear methods by more than one order of magnitude (compare the red dots with magenta triangles). At higher $\Delta t$, all methods seem to have a slope close to one. However, at small $\Delta t$, the methods with linear interpolation plateau as they get dominated by interpolation errors.

Considering that GBLA is an iterative algorithm, different sources of error are intertwined, and tracking their effects separately is complicated. However, from the numerical investigation of global errors in figure 5, we can draw several insightful conclusions. Different variations of GBLA converge robustly by lowering $\Delta x$ and $\Delta t$. The GBLA method is more affected by the interpolation error than the particle tracking, which can be mitigated by using cubic interpolation (or other higher-order schemes). If in some applications linear interpolations are selected to lighten the computations, the choice of (2.8) over (2.7) substantially improves the result.

3.2.2. Conservation of potential vorticity in shallow water

To put our method to a more challenging test, we consider a fully developed turbulent flow evolving under the rotating shallow-water equations on flat topography,

(3.6a)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{f} \times \boldsymbol{u} ={-}g \boldsymbol{\nabla} h, \end{gather}$$

(3.6b)$$\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (h \boldsymbol{u}) = 0, \end{gather}$$

where $h$ the height field, $g$ the gravitational acceleration and $\boldsymbol {f} = f \hat {\boldsymbol {z}}$ the Coriolis parameter. The existence of vortices of different sizes and sharp fronts in turbulent shallow water systems allow us to better examine the strengths and weaknesses of our method. For this flow, the PV,

(3.7)\begin{equation} q = \frac{\boldsymbol{\nabla} \times \boldsymbol{u} + f}{h}, \end{equation}

is conserved on each particle trajectory. This material conservation can be used in a similar way to the previous section to assess the accuracy of Lagrangian averaging; the LM PV expressed in terms of final position should be equal to the instantaneous PV at the end of the averaging interval. Selecting $f = 0.1$, $g = 0.1$ and the mean height $H = 1$, we solve the nonlinear shallow water equation using a similar pseudo-spectral method with the same computational resolution, domain size and boundary condition as in § 3.2.1. A viscous damping of the form $\nu \nabla ^2$ with $\nu = 8 \times 10^{-7}$ is included for numerical stability. We initialise the simulations with a fully developed turbulent flow which is initially in a geostrophic balance (see the initial PV in figure 6b). The time step of the simulation is $\delta t = 5 \times 10^3$, and that of GBLA $\Delta t = 2 \times 10^3$. The instantaneous PV at $t=120$ in figure 6(a) is remarkably close to the LM PV over the period of $[0,\ 120]$ in figure 6(c), where small-scale features are also captured by GBLA (see figure 6d for the difference between the two). The small disagreement is partly caused by the viscous damping that makes the conservation of PV less accurate. This is particularly the case for the sharp fronts in figure 6(c) which are dampened in figure 6(a).

Figure 6.

(a) The potential vorticity (PV) at the end of the averaging interval $t=120$. (b) The PV at $t=0$. (c) The LM PV from $t=0$ to $t=120$ using GBLA with a cubic interpolation scheme. (d) Absolute value of the difference between (a,c). Panels (a–c) share the same colour bar, which is different from that of panel (d).

3.3. Stokes terms for standing wave in shallow water

The difference between the Lagrangian and Eulerian mean is termed Stokes correction (see e.g. Andrews & McIntyre Reference Andrews and McIntyre1978; Bühler Reference Bühler2009) and can be written as

(3.8)\begin{equation} \overline{\phi}^{ {\hspace{0.5mm}\small S}} = \overline{\phi}^{{\hspace{0.5mm}\small L}} - \bar{\phi} = \overline{(\boldsymbol{{\xi}} \boldsymbol{\cdot} \boldsymbol{\nabla}) \phi'} + O(\epsilon^2), \quad \phi' = \phi(\boldsymbol{x},t) - \bar{\phi}, \end{equation}

where the first term in the right-hand side is the leading-order approximation for this quantity. Here $\epsilon$ is the amplitude of the wave term $\phi '$ after proper scaling, where we assumed that $\epsilon \ll 1$. We consider a monochromatic standing wave of the form

(3.9)\begin{align} h'(x,y,t) = A \cos(\omega t) \hat{h}(x,y), \quad \hat{h}(x,y) = \cos(kx) \cos(ly), \end{align}

(3.10)\begin{align} \boldsymbol{u}'(x,y,t) = \frac{g}{\omega^2 - f^2}( - \omega \sin(\omega t) \boldsymbol{\nabla} \hat{h} + \cos(\omega t) \boldsymbol{f} \times \boldsymbol{\nabla} \hat{h}), \end{align}

where $A$ is a small amplitude, is a solution of linear rotating shallow water equations provided that $\omega$ satisfies the dispersion relation

(3.11)\begin{equation} \omega = \sqrt{f^2+g H (l^2 + k^2)}. \end{equation}

Thomas, Bühler & Smith (Reference Thomas, Bühler and Smith2018) showed this standing wave induces a mean flow characterised by the Stokes terms

(3.12a,b)\begin{align} \overline{h}^{ {\hspace{0.5mm}\small S}} = \frac{1}{2} \frac{g}{\omega^2 - f^2} \left| \boldsymbol{\nabla} \hat{h} \right|^2, \quad \overline{\boldsymbol{u}}^{ {\hspace{0.5mm}\small S}} ={-}\left( \frac{g}{\omega^2 - f^2} \right)^2 \boldsymbol{f} \times \boldsymbol{\nabla} \left( \frac{1}{2} \left| \boldsymbol{\nabla} \hat{h} \right|^2 + \frac{\omega^2 - f^2}{4 g H} \hat{h}^2 \right). \end{align}

The above analytical expressions provide a good test case for GBLA. We perform a numerical simulation for

(3.13a–c)\begin{equation} A = 0.001, \quad f=g=H=1, \quad k = l =1, \end{equation}

and compare $\overline {h}^{{\hspace {0.5mm}\small L}} -\bar {h}$ and $\overline {\boldsymbol {u}}^{{\hspace {0.5mm}\small L}} -\bar {\boldsymbol {u}}$ calculated by the grid-based averaging method with the analytical expressions in (3.12a,b). There is an excellent agreement between these two sets of quantities in figure 7.

Figure 7.

Stokes terms computed by GBLA (a–c) and from analytical expression in (3.12a,b) (d–f): (a,d) $\overline {u}^{ {\hspace {0.5mm}\small S}}$; (b,e) $\overline {v}^{ {\hspace {0.5mm}\small S}}$; (c, f) $\overline {h}^{ {\hspace {0.5mm}\small S}}$. The averaging interval is set to $2{\rm \pi} /\omega$.