1. Introduction
The motions of the ocean and atmosphere involve processes with a wide range of spatial and temporal scales. In particular, fast internal waves propagate throughout these stratified, rotating fluids, interacting with slower eddies and currents. Internal waves play a key role in the momentum and energy budgets of the ocean and atmosphere, forcing the mean flow, transferring energy between large and small scales and causing turbulent mixing when they break (Naveira Garabato et al. Reference Garabato, Alberto, Nurser, Scott and Goff2013; Waterhouse Reference Waterhouse2014; MacKinnon Reference MacKinnon2017; Whalen, MacKinnon & Talley Reference Whalen, MacKinnon and Talley2018; Shakespeare & Hogg Reference Shakespeare and Hogg2019; Whalen et al. Reference Whalen, de, Casimir, Alberto, Klymak, MacKinnon and Sheen2020). Understanding and quantifying the role of internal waves in shaping the larger-scale circulation is key to designing accurate climate models, since their spatial scales cannot be directly resolved and must instead be parameterised. Wave-resolving numerical simulations are an important tool for understanding the physics that underpins these parameterisations.
A primary challenge in the study of interactions of internal waves with the non-wave flow in numerical wave-resolving simulations is separating the processes so that they can be quantified and their physics understood. We hereafter refer to the non-wave flow as the mean flow, with the understanding that this does not place any restriction on the temporal or spatial scales of the mean flow at this stage. Many different methods have been proposed for separating waves from the mean flow, often relying on defining a ‘balanced’ mean flow based on dynamical considerations (for example, geostrophic balance in the limit of small Rossby number, or its higher-order variants (Vanneste Reference Vanneste2013; Vallis Reference Vallis2017)), and taking the wave component of the flow to be the ‘unbalanced’ residual (Bühler Reference Bühler2014). However, recent studies have highlighted the importance of mean flow regimes that are
$O(1)$
in the Rossby number, indicating that inertial and rotational forces on the flow are comparable, and may therefore be ‘unbalanced’ (McWilliams Reference McWilliams2016; Taylor & Thompson Reference Taylor and Thompson2023). These flows are termed submesoscale in the ocean, and mesoscale in the atmosphere.
Submesoscale currents and internal waves are both typically energetic in the important surface and bottom boundary layers of the ocean, and despite our inability to directly capture their effects in climate models, they are starting to be regularly resolved in realistic, regional high-resolution numerical models (e.g. Nagai et al. Reference Nagai, Tandon, Kunze and Mahadevan2015; Bachman et al. Reference Bachman, Taylor, Adams and Hosegood2017; Su et al. Reference Su, Wang, Klein, Thompson and Menemenlis2018; Baker et al. Reference Baker, Mashayek, Garabato and Alberto2023). Understanding the physics of their interactions (to ultimately inform parameterisations in coarser models) has therefore become a topic of significant recent interest (e.g. Tedesco et al. Reference Tedesco, Baker, Naveira Garabato, Mazloff, Gille, Caulfield and Mashayek2023; Barkan, Srinivasan & McWilliams Reference Barkan, Srinivasan and McWilliams2024; Thomas et al. Reference Thomas, Moum, Qu, Hilditch, Kunze, Rainville and Lee2024), requiring an effective way to separate the wave-like part of the flow from other motions. In particular, to study wave generation, propagation and mixing it is important to know the wave-like part of the flow, not only the mean flow. This task is more challenging in practice than only finding the mean flow, since the waves are often lower amplitude than the mean flow and therefore more easily contaminated by imperfect decompositions.
Averaging techniques based on spatial or temporal scales are often used to define a mean state, with the wave component of the flow defined as the perturbation from this mean. In particular, weighted averages can be used to control the temporal frequencies or spatial wavenumbers that constitute the mean state (see § 2.1). We refer to these weighted averages as filters, but keep in mind that filtering is just a special case of averaging.
Filtering on wavenumber or frequency also presents problems for separating waves from the mean flow. Internal waves can have wavelengths ranging from hundreds of metres to hundreds of kilometres, often overlapping in spatial scales with motions such as oceanic submesoscales. Moreover, whilst internal waves are often considered ‘fast’ and the non-wave flow ‘slow’, this is a simplification. Internal waves in geophysical flows have an intrinsic frequency greater than
$f$
, where
$f$
is the Coriolis parameter and quantifies the rate of the Earth’s rotation (although this can be modified to an ‘effective’ Coriolis parameter by background vorticity or baroclinicity of the flow; Kunze Reference Kunze1985; Whitt & Thomas Reference Whitt and Thomas2013). However, the intrinsic frequency is the frequency in the frame of the flow, rather than the rest frame. Due to Doppler shifting of internal waves by the mean flow, there is no such frequency constraint on internal waves in the rest frame. Indeed, an important class of internal waves in the ocean and atmosphere is steady, topographically generated lee waves, which have zero frequency in the fluid’s rest frame.
Furthermore, when the waves have non-negligible amplitude they displace the mean flow with the wave frequency, so that the non-wave component of a flow found by temporally filtering in the rest frame is ‘blurred’ by the presence of waves (Kafiabad & Vanneste 2023; hereafter KV23). These two effects often render temporal filtering in the rest frame (Eulerian temporal filtering) problematic.
A proposed solution to these filtering difficulties is to perform the temporal filter not in the rest frame, but in the frame moving with the flow. This has been termed Lagrangian filtering (Nagai et al. Reference Nagai, Tandon, Kunze and Mahadevan2015; Shakespeare & Hogg Reference Shakespeare and Hogg2017). The key assumption in the geophysical context is that internal waves can be defined by a super-inertial (
${\gt}f$
) intrinsic frequency (Polzin & Lvov Reference Polzin and Lvov2011), although
$f$
can be replaced by an effective inertial frequency when the background vorticity is strong (Kunze Reference Kunze1985; Rama, Shakespeare & Hogg Reference Rama, Shakespeare and Hogg2022). Whilst Lagrangian filtering in numerical simulations is a relatively recent development, the power of Lagrangian mean theories has been well known since the pioneering theoretical work of Bretherton (Reference Bretherton1971), Soward (Reference Soward1972) and Andrews & McIntyre (Reference Andrews and McIntyre1978). In particular, the development of generalised Lagrangian mean (GLM) theory by Andrews & McIntyre (Reference Andrews and McIntyre1978) showed that a tractable definition of Lagrangian means is available and, when applied to the equations of fluid motion, leads to simple and natural equations for the Lagrangian mean that are not available for the Eulerian mean (Andrews & McIntyre Reference Andrews and McIntyre1978; Bühler Reference Bühler2014; Gilbert & Vanneste Reference Gilbert and Vanneste2018; Kafiabad, Vanneste & Young Reference Kafiabad, Vanneste and Young2021; Gilbert & Vanneste Reference Gilbert and Vanneste2025).
Despite the many benefits of Lagrangian temporal averaging (and, as a special case, Lagrangian filtering) over Eulerian averaging, it is not typically used in the analysis of numerical simulations due to computational challenges: numerical simulations are usually Eulerian in nature, with data being defined at fixed spatial grid points. Previous approaches have used particle tracking methods, whereby synthetic passive particles are seeded in a numerical simulation and advected by the flow velocities, either during the simulation or afterwards using saved data. The scalar fields recorded at the particle locations can then be used to formulate Lagrangian averages (Nagai et al. Reference Nagai, Tandon, Kunze and Mahadevan2015; Shakespeare & Hogg Reference Shakespeare and Hogg2017, Reference Shakespeare and Hogg2018, Reference Shakespeare and Hogg2019; Bachman et al. Reference Bachman, Shakespeare, Kleypas, Castruccio and Curchitser2020). Various difficulties with this approach include the computational expense of advecting particles and the potential of particles to cluster in certain parts of the domain, making a domain-wide Lagrangian average subject to potentially inaccurate interpolations. A recent open-source Python package developed by Shakespeare et al. (Reference Shakespeare, Gibson, Hogg, Bachman, Keating and Velzeboer2021) overcomes the latter of these difficulties by performing the particle tracking on offline simulation data and seeding particles at the midpoint of the interval of interest, from which they are tracked back and forth. This method has been successfully used to filter internal waves in a number of studies, but relies on large amounts of high spatial and temporal resolution simulation data being saved and processed (e.g. Shakespeare et al. Reference Shakespeare, Gibson, Hogg, Bachman, Keating and Velzeboer2021; Baker & Mashayek Reference Baker and Mashayek2022; Tedesco et al. Reference Tedesco, Baker, Naveira Garabato, Mazloff, Gille, Caulfield and Mashayek2023; Jones et al. Reference Jones, Xiao, Abernathey and Smith2023).
An alternative method for finding the Lagrangian mean was recently proposed by KV23, following on from a previous grid-based method for the same procedure (Kafiabad Reference Kafiabad2022). They showed that it is possible to define partial Lagrangian mean fields that lead to a set of partial differential equations (PDEs) that only depend on the current simulation time. These PDEs can be solved alongside the governing equations of the flow over the averaging interval, after which the final value of the partial Lagrangian means is equal to the full Lagrangian mean of interest. For each interval over which the Lagrangian mean equations are solved, only one instance of the full Lagrangian mean is computed, so the Lagrangian mean can either be found at a coarse temporal resolution, or multiple sets of the PDEs can be solved simultaneously to achieve a higher temporal resolution. This method allows the Lagrangian mean to be found on-the-fly, with no expensive data writing, storage or post-processing required. The Lagrangian mean equations can be solved with the same scheme as the governing equations of the flow.
The Lagrangian mean found using the KV23 approach is the special case of an unweighted mean over a finite interval, often referred to as a ‘top-hat’ mean. However, in order to control the frequencies that are retained by the mean field, i.e. to apply a Lagrangian frequency filter, a weighted mean is needed. In this work, we extend the method of KV23 to one for a general convolutional weighted mean. This allows us to perform Lagrangian filtering on-the-fly in a numerical simulation without particle tracking. We present three strategies for this purpose – two of which were previously presented by KV23, and a third, new strategy that avoids some issues associated with the other two.
Once the Lagrangian mean fields are computed, it is possible to define the corresponding perturbation fields in a number of ways. We therefore also carefully present several different Lagrangian wave–mean decompositions and their properties. Whilst the motivation here is to enable identification of internal inertia–gravity waves, this method is flexible in that any intrinsic frequency criterion can be used to define the wave-like perturbations. Although we focus on geophysical flows, our method can be used for any multi-time-scale flows, such as those in astrophysical or biological fluids.
This paper is structured as follows. In § 2, we introduce the weighted Lagrangian mean, some of its important properties and consider desirable forms of the weight function. In § 3, we derive the on-the-fly method for solving for the Lagrangian mean. In § 4 we introduce a rotating shallow water model that we use as a test bed for the Lagrangian mean computation, and in § 5 we show results of solving the Lagrangian mean equations alongside this model for the various strategies. Then, in § 6 we return to a more theoretical look at how a wave–mean decomposition should be defined, before presenting results of these wave–mean decompositions in the shallow water model in § 7. We discuss potential errors in § 8, and our methods and results in § 9.
2. Lagrangian mean formulation
2.1. Weighted averages and frequency filters
We begin by defining a standard weighted time average at the reference time
$t^{\ast}$
of some scalar quantity
$g(t)$
as
\begin{equation} \bar {g}(t^{\ast}) = \int _{-\infty }^{\infty } g(s)F(s,t^{\ast}) \,\mathrm {d} s, \end{equation}
where
$F(s,t^{\ast})$
is some weight function. We can design this weight function to act as a frequency filter on g at time
$t^{\ast}$
. Consider the Fourier transform of g, given by
\begin{equation} \hat {g}(\omega ) = \int _{-\infty }^{\infty }g(s)\mathrm {e}^{- \mathrm {i} \omega s}\,\mathrm {d} s, \end{equation}
then the frequency filtered scalar
$\bar {g}$
at time
$t^{\ast}$
is given by
\begin{equation} \bar {g}(t^{\ast}) = \frac {1}{2\unicode{x03C0} }\int _{-\infty }^{\infty } \hat {G}(\omega )\hat {g}(\omega )\mathrm {e}^{ \mathrm {i} \omega t^{\ast}} \, \mathrm {d}\omega , \end{equation}
where
$\hat {G}(\omega )$
weights certain frequencies. For example, when
$\hat {G}(\omega ) = 1$
in
$[-\omega _c,\omega _c]$
and is zero otherwise,
$\bar {g}(t^{\ast})$
is low-pass filtered with cutoff frequency
$\omega _c$
.
Writing (2.3) in the form of (2.1) gives
\begin{equation} \bar {g}(t^{\ast}) =\int _{-\infty }^{\infty }g(s)G(t^{\ast}-s)\,\mathrm {d} s, \end{equation}
where
\begin{equation} G(t) = \frac {1}{2\unicode{x03C0} }\int _{-\infty }^\infty \hat {G}(\omega ) \mathrm {e}^{\mathrm {i} \omega t}\, \mathrm {d} \omega , \end{equation}
thus using the convolutional weight function
$F(s,t^{\ast}) = G(t^{\ast} - s)$
in (2.1) gives the frequency filter of g at time
$t^{\ast}$
. The function
$G(t)$
can also be described as an impulse response, and
$\hat {G}(\omega )$
as the corresponding frequency response of the filter. If
$G_{\mathrm {LP}}$
describes a low-pass filter with cutoff
$\omega _c$
described above, then
\begin{equation} G_{\mathrm {LP}}(t) = \frac {\sin (\omega _c t)}{\unicode{x03C0} t}, \end{equation}
and the top-hat mean over an interval of length
$2T$
is given by
\begin{equation} G_{\mathrm {TH}}(t) = (H(t+T) - H(t-T))/2T, \end{equation}
where
$H(\cdot )$
is the Heaviside step function. We hereafter consider convolutional weight functions of the form
$F(s,t^{\ast}) = G(t^{\ast} - s)$
.
2.2. Lagrangian averaging
If a time average is calculated at a fixed point in space, then it is an Eulerian time average. If it is instead calculated along the trajectory of a particle travelling with the fluid velocity
$\boldsymbol {u}(\boldsymbol {x},t)$
, it is a Lagrangian time average.
We define a flow map
$\boldsymbol {\varphi }(\boldsymbol {a},t)$
, which gives the position of a particle labelled by
$\boldsymbol {a}$
at time
$t$
. The label
$\boldsymbol {a}$
could be taken to be the position of the particle at time
$t=0$
, so that
$\boldsymbol {\varphi }(\boldsymbol {a},0) = \boldsymbol {a}$
, and in general, we think of
$\boldsymbol {\varphi }, \boldsymbol {a} \in \mathbb {R}^2$
or
$\mathbb {R}^3$
. Following KV23, we then define the weighted Lagrangian mean flow map as
\begin{equation} \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a} ,t^{\ast}) = \int _{-\infty }^{\infty } G(t^{\ast}-s)\boldsymbol {\varphi }(\boldsymbol {a},s)\,\mathrm {d} s, \end{equation}
where
$t^{\ast}$
is the time to which the Lagrangian mean is assigned. Following Gilbert & Vanneste (Reference Gilbert and Vanneste2025) we use the double-bar notation to avoid confusion with the straightforward Eulerian average. We then define the weighted generalised Lagrangian mean of a scalar field
$f(\boldsymbol {x},t)$
by
\begin{equation} \overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast}) = \int _{-\infty }^{\infty } G(t^{\ast}-s)f(\boldsymbol {\varphi }(\boldsymbol {a},s),s) \, \mathrm {d} s. \end{equation}
For comparison, we also define the corresponding Eulerian mean
\begin{equation} \overline {f}^{E}(\boldsymbol {x},t^{\ast}) = \int _{-\infty }^{\infty } G(t^{\ast}-s)f(\boldsymbol {x},s) \, \mathrm {d} s. \end{equation}
2.3. Defining a valid weight function
We require that the weight function
$G(t^{\ast}-s)$
satisfies the normalisation condition
\begin{equation} \int _{-\infty }^{\infty } G(t^{\ast}-s)\,\mathrm {d} s = 1. \end{equation}
From (2.8), this is equivalent to requiring that the mean position of a stationary particle (i.e. when flow velocity
$\boldsymbol {u} = 0$
so that the flow map
$\boldsymbol {\varphi }$
is independent of time) is its own position. Without this natural assumption, our methods for finding
$\overline {f}^{L}$
in a periodic domain break down (see discussion of suitable domains and boundary conditions in § 3.6).
Writing
$G$
in terms of its Fourier transform
$\hat {G}$
(defined in (2.5)) shows that (2.11) is equivalent to
$\hat {G}(0) = 1$
. The filter must therefore include the zero frequency. The filters described in (2.6) and (2.7) satisfy this criterion, but a high-pass filter, for example, could not be used to define a mean flow. However, we could define a weight function that removes some frequencies
$\omega$
in a specified interval
$0 \lt |\omega _1| \lt |\omega | \lt |\omega _2|$
by setting
\begin{equation} \hat {G}(\omega ) = \begin{cases} 0 , \hspace {1cm} |\omega _1| \lt |\omega | \lt |\omega _2|, \\ 1 , \hspace {1cm} \mathrm {otherwise}. \end{cases} \end{equation}
We later show an example of this filter in figure 7.
We would also like the weight function to be such that
$\overline {f}^{L}$
(defined in (2.9)) satisfies a property that we expect of a mean, namely, that the mean is unchanged by reapplying the averaging operation
\begin{equation} \overline{{\,\overline {f}}^{L}}^{L} (\boldsymbol {x},t^{\ast})\, \overset {?}{=} \overline {f}^{L} (\boldsymbol {x},t^{\ast}). \end{equation}
Since the Lagrangian mean of a scalar depends on the flow with which it is advected, there is some ambiguity in the notation on the left-hand side of (2.13). We formalise this statement in Appendix A, and note here that the Lagrangian mean of
$\overline {f}^{L}$
is taken with respect to the Lagrangian mean flow defined by the map
$\bar {\bar {\boldsymbol {\varphi }}}$
.
Equation (2.13) is satisfied only when
$\hat {G}(\omega ) = 0$
or
$1$
, or some piece-wise combination of each. A proof of this is given in Appendix A. This is equivalent to requiring that
$\hat {G}$
represents a perfect band-pass filter. In this case, the Lagrangian filtered field
$\overline {f}^{L}$
behaves as we hope a ‘mean’ field should, in that it contains no (Lagrangian) high frequencies. In practice, filters rarely exactly satisfy the condition (2.13). Perfect band-pass filters tend to suffer from spectral ringing at the cutoff frequency, so other imperfect filters such as the Butterworth filter are often used instead (Rama et al. Reference Rama, Shakespeare and Hogg2022). Here, the wave frequency considered in the numerical model in § 4 is not close to the cutoff frequency, so ringing is not an issue. We therefore only consider perfect band-pass filters so that the mean field is not expected to contain any wave signal. Using a fourth-order Butterworth filter gives indistinguishable results in our case, although it does allow shorter averaging intervals (see discussion of figure 4), which makes it preferable in practice, even though (2.13) is not perfectly satisfied.
We note here that there is no assumption of time-scale separation between the ‘slow’ mean flow and the ‘fast’ motions to be filtered. The original formulation of GLM theory by Andrews & McIntyre (Reference Andrews and McIntyre1978) defines the Lagrangian average in an abstract way to apply to ensemble averages. To apply it to temporal averages such as the one we compute requires an assumption that the mean flow is ‘frozen’ during the averaging operation (Bühler Reference Bühler2014), and this is explained further with an example illustrating the difference between the formulations in Appendix B.
2.4. Lagrangian mean velocity
By definition, the flow map
$\boldsymbol {\varphi }$
satisfies
\begin{equation} \boldsymbol {u}(\boldsymbol {\varphi }(\boldsymbol {a},t),t) = \frac {\partial \boldsymbol {\varphi }}{\partial t}(\boldsymbol {a},t), \end{equation}
where
$\boldsymbol {u}$
is the fluid velocity. The Lagrangian mean velocity
$\bar {\bar {\boldsymbol {u}}}$
is then defined to be the velocity of a particle moving along a Lagrangian mean trajectory, that is,
\begin{equation} \bar {\bar {\boldsymbol {u}}}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast}) = \frac {\partial \bar {\bar {\boldsymbol {\varphi }}}}{\partial t^{\ast}}(\boldsymbol {a},t^{\ast}). \end{equation}
However, another velocity
$\overline {\boldsymbol {u}}^{L}$
can be defined by taking the Lagrangian mean of each component of the velocity
$\boldsymbol {u}$
treated as scalars (see Gilbert & Vanneste (Reference Gilbert and Vanneste2018) and Gilbert & Vanneste (Reference Gilbert and Vanneste2025) for other, more geometric definitions of
$\overline {\boldsymbol {u}}^{L}$
). The averaging operation is such that
$\bar {\bar {\boldsymbol {u}}} = \overline {\boldsymbol {u}}^{L}$
for the class of convolutional weight functions considered here. This is a special case of the more general result
\begin{equation} \bar {\bar {D}}\overline {f}^{L} = \overline {\,D f\,}^{L}, \end{equation}
where
\begin{align} D &\equiv \frac {\partial }{\partial t} + \boldsymbol {u}\boldsymbol {\cdot } \nabla , \end{align}
\begin{align} \bar {\bar {D}} &\equiv \frac {\partial }{\partial t^{\ast}} + \bar {\bar {\boldsymbol {u}}}\boldsymbol {\cdot } \nabla . \end{align}
This result is shown in Appendix C, with the result
$\bar {\bar {\boldsymbol {u}}} = \overline {\boldsymbol {u}}^{L}$
found by considering
$f(\boldsymbol {x},t^{\ast}) = \boldsymbol {x}$
. Equation (2.16) is one of the most powerful results of the Lagrangian formalism – it means that material conservation laws and scalar transport relations are inherited naturally by the corresponding Lagrangian means (Andrews & McIntyre Reference Andrews and McIntyre1978).
We now formulate a PDE-based method for calculating these Lagrangian mean quantities, extending the method of KV23 to include a general convolutional weight function
$G(t^{\ast}-s)$
, which allows us to use a Lagrangian filter for wave–mean decomposition.
3. Formulation of on-the-fly method
KV23 developed a method for finding the top-hat Lagrangian mean
$\overline {f}^{L}$
of a scalar field
$f$
, as defined in (2.9) (with
$G(t) = (H(t+T) - H(t-T))/2T$
), by formulating equations for the ‘partial Lagrangian means’ and evolving them in a numerical simulation alongside the governing equations of the flow. Here, we re-derive this method for a weighted Lagrangian mean, although we have the specific application of a low-pass filter in mind.
KV23 presented two strategies for finding
$\overline {f}^{L}$
. Strategy 1 solves first for an auxiliary mean function, before using a remapping to recover
$\overline {f}^{L}$
, whereas strategy 2 solves directly for
$\overline {f}^{L}$
. Both of these strategies have particular advantages and disadvantages, which we discuss further later. Here, we re-derive these two strategies for a weighted mean, and present a new third strategy that circumvents some difficulties with strategies 1 and 2.
3.1. Definition of full Lagrangian means
We now approximate the average over an infinitely long interval in (2.8)–(2.9) by one over a finite interval
$[t^{\ast}-T,t^{\ast}+T]$
, which is centred on the time
$t^{\ast}$
at which the average is defined. This need not be the case, but it is the natural choice for even weight functions
$G$
, which correspond to real frequency response functions
$\hat {G}$
(see (2.5)).
The mean flow map is now defined by (cf. (2.8))
\begin{equation} \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}) = \int _{t^{\ast}-T}^{t^{\ast}+T} G(t^{\ast}-s) \boldsymbol {\varphi }(\boldsymbol {a},s)\,\mathrm {d} s, \end{equation}
where
$G(t^{\ast}-s)$
satisfies the normalisation (2.11) over the interval
$[t^{\ast}-T,t^{\ast}+T]$
.
We introduce some notation defining rearrangements of the mean scalar
$\overline {f}^{L}$
that depend on different spatial coordinates
\begin{equation} \begin{aligned}\overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast}) &= \tilde {f}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}+T),t^{\ast}) = f^{\ast}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast})\\ &= \int _{t^{\ast}-T}^{t^{\ast}+T} G(t^{\ast}-s) f(\boldsymbol {\varphi }(\boldsymbol {a},s),s)\,\mathrm {d} s . \end{aligned}\end{equation}
The variables
$\overline {f}^{L}$
,
$\tilde {f}$
and
$f^{\ast}$
all encode the Lagrangian mean of
$f$
, but they use different independent variables to do so:
$\overline {f}^{L}(\boldsymbol {x},t^{\ast})$
is the Lagrangian mean for the particle whose mean position is
$\boldsymbol {x}$
,
$\tilde {f}(\boldsymbol {x},t^{\ast})$
is the Lagrangian mean for the particle whose position at time
$t^{\ast} + T$
is
$\boldsymbol {x}$
and
$f^{\ast}(\boldsymbol {x},t^{\ast})$
is the Lagrangian mean for the particle whose position at time
$t^{\ast}$
is
$\boldsymbol {x}$
. The three functions are rearrangements of each other, that is, related via composition with (not necessarily volume-preserving) maps.
Despite simply being rearrangements of each other, each of the definitions in (3.2) will be useful to us. The variable
$\overline {f}^{L}$
is the generalised Lagrangian mean that we are looking to find, and is the true Lagrangian mean in that it satisfies properties such as (2.13) and (2.16). The variables
$\tilde {f}$
and
$f^{\ast}$
are auxiliary fields that will help us to derive
$\overline {f}^{L}$
in two of our strategies, and we will also demonstrate that
$f^{\ast}$
is useful in itself to extract the wave field. Hence, one may want to compute it in addition to
$\overline {f}^{L}$
.
We also need to define the set of maps
$\boldsymbol {\Xi }$
between each of the spatial independent variables that effect the rearrangements of
$\overline {f}^{L}$
in (3.2). For this purpose, the label ‘1’ refers to the trajectory endpoint position
$\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}+T)$
, ‘2’ refers to the trajectory mean position
$\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})$
and ‘3’ refers to the trajectory midpoint position
$\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast})$
. We use this convention because strategy ‘
$i$
’ directly finds the Lagrangian mean in (3.2) with spatial independent variable ‘
$i$
’, for
$i \in \{1,2,3\}$
. For example, strategy 1 solves directly for
$\tilde {f}$
. We define the map
$\boldsymbol {\Xi }^{i\mapsto j}$
to map from the
$i$
coordinate to the
$j$
coordinate, such that
$\boldsymbol {\Xi }^{i\mapsto j}$
is the identity map for
$i=j$
,
$ (\boldsymbol {\Xi }^{i\mapsto j} )^{-1} = \boldsymbol {\Xi }^{j\mapsto i}$
and
\begin{align} \boldsymbol {\Xi }^{1\mapsto 2}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}+T),t^{\ast}) &= \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}), \end{align}
\begin{align} \boldsymbol {\Xi }^{1\mapsto 3}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}+T),t^{\ast}) &= \boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}), \end{align}
\begin{align} \boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast}) &= \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}). \end{align}
These maps are illustrated in the schematic in figure 1.
Figure 1. Schematic of a particle trajectory (black) with label
$\boldsymbol {a}$
in the interval
$[t^{\ast}-T,t^{\ast} + T]$
, with positions labelled by the flow map
$\boldsymbol {\varphi }(\boldsymbol {a},t)$
. The mean particle trajectory on the same interval is shown in blue, with positions labelled by the mean flow map
$\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})$
. Red arrows indicate the maps
$\boldsymbol {\Xi }^{i \mapsto j}$
from position
$i$
to position
$j$
, where position 1 is the trajectory endpoint
$\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast} +T)$
, position 2 is the trajectory mean
$\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})$
and position 3 is the trajectory midpoint
$\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast})$
.
3.2. Definition of partial Lagrangian means
As in KV23, we now define a corresponding set of ‘partial’ Lagrangian mean fields, that is, fields obtained by carrying out the averaging integration from
$t^{\ast} - T$
to some
$t \lt t^{\ast} + T$
. By finding PDEs for these partial fields and evolving them over the averaging interval alongside the governing equations for the flow, the full Lagrangian mean fields in § 3.1 are obtained when
$t = t^{\ast} + T$
. The subscript
$p$
always denotes a ‘partial’ field, and these fields evolve with time
$t$
, while the time
$t^{\ast}$
at which the Lagrangian mean is assigned is a fixed parameter. In the definitions of the partial fields, we drop the dependence on
$t^{\ast}$
for readability, since everything in this section refers to one fixed averaging time
$t^{\ast}$
.
First, we define a partial mean flow map to correspond to (3.1), namely
\begin{equation} \bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t) = \int _{t^{\ast}-T}^t G(t^{\ast}-s) \boldsymbol {\varphi }(\boldsymbol {a},s)\,\mathrm {d} s + \boldsymbol {\varphi }(\boldsymbol {a},t)\left (1 - \int _{t^{\ast}-T}^t G(t^{\ast}-s)\,\mathrm {d} s\right ), \end{equation}
so that
$\bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t^{\ast}+T) = \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})$
. This particular form of
$\bar {\bar {\boldsymbol {\varphi }}}_p$
(in particular the second term which vanishes when
$t = t^{\ast} + T$
) is needed for a similar reason to that discussed in § 2.3, namely that the partial mean position of a stationary particle should be the position itself, so that the image of the partial mean flow map is the same as that of the flow map itself.
We then define the partial equivalents of (3.2)
\begin{equation} \overline {f}^{L}_p(\bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t),t) = \tilde {f}_p(\boldsymbol {\varphi }(\boldsymbol {a},t),t) = f^{\ast}_p({\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t),t) = \int _{t^{\ast}-T}^t G(t^{\ast}-s) f(\boldsymbol {\varphi }(\boldsymbol {a},s),s)\,\mathrm {d} s, \end{equation}
where a second term corresponding to that used in (3.6) is not necessary here and is omitted for convenience. The new partial coordinate corresponding to
$\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast})$
is
${\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t)$
, given by
\begin{equation} {\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t) = \begin{cases} \boldsymbol {\varphi }(\boldsymbol {a},t), & t \lt t^{\ast}, \\ \boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}), & t \geqslant t^{\ast}. \end{cases} \end{equation}
This form of
${\boldsymbol {\varphi }}^{\ast}_p$
is necessary because using
$\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast})$
as the coordinate from the beginning would violate causality. Other definitions, such as
${\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t) = \boldsymbol {\varphi }(\boldsymbol {a},t/2)$
, although resulting in the desired full mean definition
$f^{\ast}$
(see (3.2)), would not allow the subsequent evolution equations to depend only on fields at the current time.
The definitions in (3.7) ensure that the full Lagrangian means defined in (3.2) are recovered by setting
$t = t^{\ast} + T$
in (3.7).
We also define partial mean equivalents of (3.3)–(3.5) (where as before,
$ (\boldsymbol {\Xi }^{i\mapsto j}_p )^{-1} = \boldsymbol {\Xi }^{j\mapsto i}_p$
)
\begin{align} \boldsymbol {\Xi }^{1\mapsto 2}_p(\boldsymbol {\varphi }(\boldsymbol {a},t),t) &= \bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t), \end{align}
\begin{align} \boldsymbol {\Xi }^{1\mapsto 3}_p(\boldsymbol {\varphi }(\boldsymbol {a},t),t) &= {\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t), \end{align}
\begin{align} \boldsymbol {\Xi }^{2\mapsto 3}_p(\bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t),t) &= {\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t). \end{align}
We now have all of the notation necessary to derive three separate strategies for finding
$\overline {f}^{L}$
and/or
$f^{\ast}$
, which may be sufficient (see § 6). We briefly summarise each strategy here. In table 1, we summarise the dependent fields of each PDE to be solved for the three strategies.
Table 1. Fields to be solved for in each of the three strategies presented.
-
(i) Strategy 1: Solve for
$\tilde {f} (\boldsymbol {x},t^{\ast})$
and the map
$\boldsymbol {\Xi }^{1\mapsto 2}(\boldsymbol {x},t^{\ast})$
, then find
$\overline {f}^{L}(\boldsymbol {x},t^{\ast}) = \tilde {f}((\boldsymbol {\Xi }^{1\mapsto 2})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast})$
by a final interpolation step. The variable
$\ f^{\ast}(\boldsymbol {x},t^{\ast}) = \tilde {f}((\boldsymbol {\Xi }^{1\mapsto 3})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast})$
can also be found by solving for
$\boldsymbol {\Xi }^{1\mapsto 3}(\boldsymbol {x},t^{\ast})$
. -
(ii) Strategy 2: Solve directly for
$\overline {f}^{L} (\boldsymbol {x},t^{\ast})$
, which requires also solving for the map
$\boldsymbol {\Xi }^{2\mapsto 1}(\boldsymbol {x},t^{\ast})$
. The variable
$f^{\ast}(\boldsymbol {x},t^{\ast}) = \overline {f}^{L} ((\boldsymbol {\Xi }^{2\mapsto 3})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast})$
can also be found by solving for
$\boldsymbol {\Xi }^{2\mapsto 3}(\boldsymbol {x},t^{\ast})$
. -
(iii) Strategy 3: Solve directly for
$ f^{\ast} (\boldsymbol {x},t^{\ast})$
, which requires also solving for the map
$\boldsymbol {\Xi }^{3\mapsto 1}(\boldsymbol {x},t^{\ast})$
. The variable
$\overline {f}^{L}(\boldsymbol {x},t^{\ast}) = f^{\ast} ((\boldsymbol {\Xi }^{3\mapsto 2})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast})$
can be found by solving for
$\boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {x},t^{\ast})$
.
3.3.
Strategy 1: solve at partial trajectory endpoint
$\boldsymbol {\varphi }(\boldsymbol {a},t)$
The first strategy consists of solving for
$\tilde {f}(\boldsymbol {x},t^{\ast})$
, the Lagrangian mean along the trajectory of a particle whose position
$\boldsymbol {x} = \boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}+T)$
is at the trajectory endpoint (as defined in (3.2)).
Taking the time derivative of (3.7) at fixed
$\boldsymbol {a}$
, and using the dummy variable
$\boldsymbol {x} = \boldsymbol {\varphi }(\boldsymbol {a},t)$
, we find
\begin{equation} \frac {\partial \tilde {f}_p}{\partial t}(\boldsymbol {x},t) + \boldsymbol {u}(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla \tilde {f}_p(\boldsymbol {x},t) = f(\boldsymbol {x},t)G(t^{\ast} - t), \end{equation}
where we have used
$\boldsymbol {u}(\boldsymbol {\varphi }(\boldsymbol {a},t),t) = {\partial \boldsymbol {\varphi }}/{\partial t}(\boldsymbol {a},t)$
by definition of the flow map.
Equation (3.12), along with this initial condition
$\tilde{f}_p(\mathbf{x},t^{\ast}-T) = 0$
(see (3.7)), can be solved alongside the governing equations to find
$\tilde {f}(\boldsymbol {x},t^{\ast}) = \tilde {f}_p(\boldsymbol {x},t^{\ast}+T)$
. However, we then need to map to the
$\bar {\bar {\boldsymbol {\varphi }}}$
coordinates to find
$\overline {f}^{L}$
. We therefore differentiate the definition of
$\boldsymbol {\Xi }^{1\mapsto 2}_p$
in (3.9), and use the definition (3.6) of
$\bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t)$
(setting
$\boldsymbol {x} = \boldsymbol {\varphi }(\boldsymbol {a},t)$
as before) to give
\begin{equation} \frac {\partial \boldsymbol {\xi }^{1\mapsto 2}_p}{\partial t}(\boldsymbol {x},t) + \boldsymbol {u}(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla \boldsymbol {\xi }^{1\mapsto 2}_p (\boldsymbol {x},t) = - \boldsymbol {u}(\boldsymbol {x},t)\int _{t^{\ast}-T}^t G(t^{\ast}-s)\,\mathrm {d} s, \end{equation}
where the perturbation map
$\boldsymbol {\xi }^{1\mapsto 2}_p$
is defined by
\begin{equation} \boldsymbol {\xi }^{1\mapsto 2}_p(\boldsymbol {x},t) = \boldsymbol {\Xi }^{1\mapsto 2}_p(\boldsymbol {x},t) - \boldsymbol {x}. \end{equation}
Initial conditions for (3.13) are given by
$\boldsymbol {\xi }^{1\mapsto 2}_p(\boldsymbol {x},t^{\ast}-T) = 0$
. Evolving (3.13) alongside (3.12) and the governing equations, we can then use
$\boldsymbol {\Xi }^{1\mapsto 2}(\boldsymbol {x},t^{\ast}) = \boldsymbol {\Xi }^{1\mapsto 2}_p(\boldsymbol {x},t^{\ast}+T) = \boldsymbol {\xi }^{1\mapsto 2}_p(\boldsymbol {x},t^{\ast}+T) + \boldsymbol {x}$
to remap
$\tilde {f}$
to
$\overline {f}^{L}$
, so that, using (3.2) and writing in terms of a dummy variable
$\boldsymbol {x}$
,
\begin{equation} \overline {f}^{L}(\boldsymbol {x},t^{\ast}) = \tilde {f}((\boldsymbol {\Xi }^{1\mapsto 2})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast}) . \end{equation}
If we would also like to know the mean field
$f^{\ast}$
defined at the flow trajectory midpoint, we can solve for
$\boldsymbol {\Xi }^{1\mapsto 3}(\boldsymbol {x},t^{\ast})$
by differentiating the definition of
$\boldsymbol {\Xi }^{1\mapsto 3}_p(\boldsymbol {x},t)$
in (3.10) and using (3.8) to find
\begin{equation} \frac {\partial \boldsymbol {\xi }^{1\mapsto 3}_p}{\partial t}(\boldsymbol {x},t) + \boldsymbol {u}(\boldsymbol {x},t) \boldsymbol {\cdot } \nabla \boldsymbol {\xi }^{1\mapsto 3}_p(\boldsymbol {x},t) = -\boldsymbol {u}(\boldsymbol {x},t)H(t-t^{\ast}), \end{equation}
where
$H(\cdot )$
is the Heaviside step function, and the perturbation map
$\boldsymbol {\xi }^{1\mapsto 3}_p$
is defined by
\begin{equation} \boldsymbol {\xi }^{1\mapsto 3}_p(\boldsymbol {x},t) = \boldsymbol {\Xi }^{1\mapsto 3}_p(\boldsymbol {x},t) - \boldsymbol {x} . \end{equation}
After (3.16) has been evolved to time
$t^{\ast}+T$
,
$f^{\ast}(\boldsymbol {x},t^{\ast})$
can then be found using
$\boldsymbol {\Xi }^{1\mapsto 3}(\boldsymbol {x},t^{\ast}) = \boldsymbol {\Xi }^{1\mapsto 3}_p(\boldsymbol {x},t^{\ast}+T) = \boldsymbol {\xi }^{1\mapsto 3}_p(\boldsymbol {x},t^{\ast}+T) + \boldsymbol {x}$
by interpolation from
\begin{equation} f^{\ast}(\boldsymbol {x},t^{\ast}) = \tilde {f}((\boldsymbol {\Xi }^{1\mapsto 3})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast}) . \end{equation}
Strategy 1 is the simplest and cheapest strategy, since the evolution equations (3.12), (3.13) and (3.16) do not involve interpolation at each time step (as will be required by strategies 2 and 3). However, in a complex flow with time scales similar to or smaller than the averaging interval, the maps
${\boldsymbol\Xi }^{1 \mapsto 2}$
and
$ {\boldsymbol\Xi }^{1 \mapsto 3}$
can be far from the identity map, making the final interpolation step inaccurate. This will be discussed further in § 8, and a case where this interpolation is too complex will be shown in figure 4. For this reason, KV23 developed a strategy 2, which avoids the need for this interpolation step.
3.4.
Strategy 2: solve at trajectory partial mean
$\bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t)$
Following KV23, in strategy 2 we solve directly for
$\overline {f}^{L}(\boldsymbol {x},t^{\ast})$
, the Lagrangian mean along the trajectory of a particle whose position
$\boldsymbol {x} = \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})$
is at the trajectory mean position (as defined in (3.2)).
Taking the derivative with respect to
$t$
of (3.6), we define
\begin{equation} \bar {\boldsymbol {u}}_p(\bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t),t) \equiv \frac {\partial \bar {\bar {\boldsymbol {\varphi }}}_p}{\partial t} = \boldsymbol {u}(\boldsymbol {\varphi }(\boldsymbol {a},t),t)\left (1 - \int _{t^{\ast}-T}^t G(t^{\ast}-s) \,\mathrm {d} s\right ). \end{equation}
We note that
$\bar {\boldsymbol {u}}_p$
is not related to the Lagrangian mean velocity
$\bar {\bar {\boldsymbol {u}}}$
defined in § 2.4. Here,
$\bar {\boldsymbol {u}}_p$
is found by taking the time derivative of
$\bar {\bar {\boldsymbol {\varphi }}}_p$
with respect to
$t$
, whereas the Lagrangian mean velocity is the derivative of
$\bar {\bar {\boldsymbol {\varphi }}}$
with respect to
$t^{\ast}$
.
Then, differentiating the definition of
$\boldsymbol {\Xi }^{2\mapsto 1}_p$
in (3.9) and letting
$\boldsymbol {x} = \bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t)$
gives
\begin{equation} \frac {\partial \boldsymbol {\xi }^{2\mapsto 1}_p}{\partial t}(\boldsymbol {x},t) + \bar {\boldsymbol {u}}_p(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla \boldsymbol {\xi }^{2\mapsto 1}_p (\boldsymbol {x},t) = \boldsymbol {u}(\boldsymbol {x} + \boldsymbol {\xi }^{2\mapsto 1}_p(\boldsymbol {x},t),t)\int _{t^{\ast}-T}^t G(t^{\ast}-s) \,\mathrm {d} s , \end{equation}
where the perturbation map
$\boldsymbol {\xi }^{2\mapsto 1}_p$
is defined by
\begin{equation} \boldsymbol {\xi }^{2\mapsto 1}_p(\boldsymbol {x},t) = \boldsymbol {\Xi }^{2\mapsto 1}_p(\boldsymbol {x},t) - \boldsymbol {x}, \end{equation}
with initial condition
$\boldsymbol {\xi }^{2\mapsto 1}_p(\boldsymbol {x},t^{\ast} - T) = 0$
, and from (3.19),
\begin{equation} \bar {\boldsymbol {u}}_p(\boldsymbol {x},t) = \boldsymbol {u}(\boldsymbol {x} + \boldsymbol {\xi }^{2\mapsto 1}_p(\boldsymbol {x},t),t)\left (1 - \int _{t^{\ast}-T}^t G(t^{\ast}-s) \,\mathrm {d} s\right ). \end{equation}
The evolution of the partial Lagrangian mean scalar
$\overline {f}^{L}_p$
can be found by taking the time derivative of (3.7), giving
\begin{equation} \frac {\partial \overline {f}^{L}_p}{\partial t}(\boldsymbol {x},t) + \bar {\boldsymbol {u}}_p(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla \overline {f}^{L}_p(\boldsymbol {x},t) = f(\boldsymbol {x} + \boldsymbol {\xi }^{2\mapsto 1}_p(\boldsymbol {x},t),t)G(t^{\ast} - t), \end{equation}
with initial condition
$\overline {f}^{L}_p(\boldsymbol {x},t^{*} - T) = 0$
.
Solving the system of equations (3.20)–(3.23) then directly gives the Lagrangian mean
$\overline {f}^{L}$
. If desired, we can also find
$f^{\ast}(\boldsymbol {x},t^{\ast})$
by solving for the map
$\boldsymbol {\Xi }^{2\mapsto 3}(\boldsymbol {x},t^{\ast})$
. Differentiating the definition of the map
$\boldsymbol {\Xi }^{2\mapsto 3}_p$
in (3.11), and setting
$\boldsymbol {x} = \bar {\bar {\boldsymbol {\varphi }}}_p(\boldsymbol {a},t)$
leads to
\begin{equation}\begin{aligned} &\frac {\partial \boldsymbol {\xi }^{2\mapsto 3}_p}{\partial t}(\boldsymbol {x},t) + \bar {\boldsymbol {u}}_p(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla \boldsymbol {\xi }^{2\mapsto 3}_p (\boldsymbol {x},t)\\ &\qquad = \boldsymbol {u}(\boldsymbol {x} + \boldsymbol {\xi }^{2\mapsto 1}_p(\boldsymbol {x},t),t)\left (\int _{t^{\ast}-T}^t G(t^{\ast}-s) \,\mathrm {d} s - H(t - t^{\ast})\right ), \end{aligned}\end{equation}
where the perturbation map
$\boldsymbol {\xi }^{2\mapsto 3}_p$
is defined by
\begin{equation} \boldsymbol {\xi }^{2\mapsto 3}_p(\boldsymbol {x},t) = \boldsymbol {\Xi }^{2\mapsto 3}_p(\boldsymbol {x},t) - \boldsymbol {x}. \end{equation}
Then,
$f^{\ast}(\boldsymbol {x},t^{\ast})$
can be found by interpolation according to
\begin{equation} f^{\ast}(\boldsymbol {x},t^{\ast}) = \overline {f}^{L}((\boldsymbol {\Xi }^{2\mapsto 3})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast}) . \end{equation}
Strategy 2 is intended to eliminate the problems with the final interpolation step in strategy 1 by solving directly at the partial Lagrangian mean position. However, this comes at the expense of a more complicated right-hand side in the evolution equations (3.23), (3.20) and (3.24), which require interpolation at each time step.
A key disadvantage of strategy 2 pertains to the boundaries of the fluid domain. Since Lagrangian variables are referenced to the trajectory mean position, the equations are posed on a moving domain that will not in general coincide with the fluid domain, making boundary conditions non-trivial – this is discussed further in § 3.6. We therefore now derive a third strategy that enables Lagrangian filtering in more complex and realistic domains.
3.5.
Strategy 3: solve at trajectory midpoint
${\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t)$
We solve directly for
$f^{\ast}(\boldsymbol {x},t^{\ast})$
, the Lagrangian mean along the trajectory of a particle whose position
$\boldsymbol {x} = \boldsymbol {\varphi }(\boldsymbol {a},t^{\ast})$
is at the trajectory midpoint (as defined in (3.2)). For this, we need to solve for the map
$\boldsymbol {\Xi }^{3\mapsto 1}$
, and also for the map
$\boldsymbol {\Xi }^{3\mapsto 2}$
if we also want to find
$\overline {f}^{L}(\boldsymbol {x},t^{\ast})$
. The derivation is similar to that for strategies 1 and 2, although the first and second halves of the interval must be considered separately. A full derivation is given in Appendix D.
Strategy 3 consists of solving (from (D1) and (D6))
\begin{equation} \frac {\partial f^{\ast}_p}{\partial t}(\boldsymbol {x},t) = G(t^{\ast} - t)f(\boldsymbol {x} + \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t),t) - H(t^{\ast}-t)\boldsymbol {u}(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla f^{\ast}_p(\boldsymbol {x},t), \end{equation}
with initial conditions
$f^{\ast}_p(\boldsymbol {x},t^{\ast} - T) = 0$
, along with (from (D3) and (D8))
\begin{equation} \frac {\partial \boldsymbol {\xi }^{3\mapsto 1}_p}{\partial t}(\boldsymbol {x},t) = H(t-t^{\ast})\boldsymbol {u}(\boldsymbol {x} + \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t),t), \end{equation}
with initial conditions
$\boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t^{\ast} - T) = 0$
. If
$\overline {f}^{L}(\boldsymbol {x},t^{\ast})$
is required, then we also solve (from (D5) and (D8))
\begin{align}\frac {\partial \boldsymbol {\xi }^{3\mapsto 2}_p}{\partial t}(\boldsymbol {x},t) &= \boldsymbol {u}(\boldsymbol {x} + \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t),t)\left (H(t-t^{\ast}) - \int _{t^{\ast}-T}^t G(t^{\ast}-s)\,\mathrm {d} s\right ) \notag \\&\quad - H(t^{\ast}-t)\boldsymbol {u}(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla \boldsymbol {\xi }^{3\mapsto 2}_p (\boldsymbol {x},t), \end{align}
with initial conditions
$\boldsymbol {\xi }^{3\mapsto 2}_p(\boldsymbol {x},t^{\ast} - T) = 0$
. Then,
$\overline {f}^{L}(\boldsymbol {x},t^{\ast})$
can be found using
$\boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {x},t^{\ast}) = \boldsymbol {\Xi }^{3\mapsto 2}_p(\boldsymbol {x},t^{\ast}+T) = \boldsymbol {\xi }^{3\mapsto 2}_p(\boldsymbol {x},t^{\ast}+T) + \boldsymbol {x}$
from
\begin{equation} \overline {f}^{L}(\boldsymbol {x},t^{\ast}) = f^{\ast}((\boldsymbol {\Xi }^{3\mapsto 2})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast}) . \end{equation}
Like strategy 1, strategy 3 requires a final interpolation step if
$\overline {f}^{L}$
is required (rather than
$f^{\ast}$
). However, this final interpolation, performed using
$\boldsymbol {\Xi }^{3\mapsto 2}$
, is likely to be much more accurate than that in strategy 1, since the trajectory mean and midpoint positions differ only by the wave perturbation (see figure 1). This will be demonstrated in figure 4.
3.6. Boundary conditions
In this study, we consider the simplest case of a doubly periodic domain. This is simple to implement as the Lagrangian mean equations for each strategy are constructed so that periodic state fields of the simulation lead to periodic Lagrangian mean fields (the normalisation condition (2.11) is essential for this to be the case). However, some of the equations can also be solved in more complex domains.
Any fluid in a domain with open (non-periodic) boundaries will contain trajectories that exit the domain, so all definitions of Lagrangian means for these trajectories will be undefined and Lagrangian mean fields cannot be calculated over the full domain. However, the equations of strategies 1 and 3 ((3.12), (3.13), (3.16) and (3.27)–(3.29)) can be straightforwardly solved in any domain with fixed boundaries. Each of the equations for Lagrangian fields in these strategies contains an advective derivative term
$\boldsymbol {u} \boldsymbol {\cdot } \nabla$
, indicating that a boundary condition is needed. However, in a fixed bounded domain with normal
$\boldsymbol {n}$
, the velocity satisfies
$\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {n} = 0$
, so the normal part of the advective term vanishes at the boundary and no boundary conditions on the Lagrangian fields are necessary. After having solved for
$\tilde {f}$
(strategy 1) or
$f^{\ast}$
(strategy 3), the final interpolation can then be carried out to find
$\overline {f}^{L}$
, although there is no guarantee that
$\overline {f}^{L}$
can be defined at every point in the domain (i.e. for a domain
$\mathcal {D}$
and some
$\boldsymbol{y} \in \mathcal {D}$
, there may not exist
$\boldsymbol {x} \in \mathcal {D}$
such that
$\boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {x},t) = \boldsymbol{y}$
).
In contrast, strategy 2 cannot easily be used in fixed bounded domains. Since the Lagrangian fields are defined on the image of the partial mean flow map
$\bar {\bar {\boldsymbol {\varphi }}}_p$
, the PDEs are posed on a domain with moving boundaries, leading in general to a free boundary problem. There is no guarantee that the Lagrangian mean position itself lies in the fluid domain (unless it is convex), or that a given location in the fluid domain is the Lagrangian mean position of some trajectory, so
$\overline {f}^{L}$
may not be defined everywhere.
Strategy 2 can, however, be straightforwardly used in domains that have at most one non-periodic dimension, along which the boundaries must align with a constant coordinate surface in that dimension (i.e. one set of straight and parallel boundaries in Euclidean space). In this case, the image of the mean flow map coincides with the fluid domain. Boundary conditions are not needed, since trajectories stay on the boundary and thus
$\boldsymbol {u}(\boldsymbol {x},t) \boldsymbol {\cdot } \boldsymbol {n} = 0 \Rightarrow \boldsymbol {u}(\boldsymbol {x} + \boldsymbol {\xi }^{2\mapsto 1}_p,t)\boldsymbol {\cdot } \boldsymbol {n} =0 \Rightarrow \boldsymbol {u}_p(\boldsymbol {x},t) \boldsymbol {\cdot } \boldsymbol {n} =0$
.
4. Numerical model
We now demonstrate our filtering approach using a single layer rotating shallow water system, which permits both geostrophic turbulence and Poincaré waves.
Similarly to KV23, we use the rotating shallow water equations in a doubly periodic domain. However, in order to have more flexibility over the chosen wavenumbers of Poincaré waves, we use the modified shallow water (MSW) equations introduced by Bühler (Reference Bühler1998). These equations were developed for the very purpose of providing a simple test bed for wave–mean decompositions, without the added complication of steepening Poincaré waves that occurs in the regular shallow water equations (Bühler Reference Bühler1998). The MSW equations behave similarly to the shallow water equations, and the equations are identical when linearised about a state of rest. In our case, we want a flow that contains a slowly varying ‘mean’ component alongside a wave field, and are agnostic to the physicality of the flow. We work with non-dimensional quantities – see KV23 for details of the non-dimensionalisation. The flow equations are
\begin{align} \frac {\partial \boldsymbol {u}}{\partial t} + \boldsymbol {u} \boldsymbol {\cdot } \nabla \boldsymbol {u} + \frac {1}{Ro}\hat {\boldsymbol {z}}\times \boldsymbol {u} &= -\frac {1}{Fr^2}\mathcal {F}(h)\nabla h, \end{align}
\begin{align} \frac {\partial h}{\partial t} + \nabla \boldsymbol {\cdot } (\boldsymbol {u} h) &= 0 , \end{align}
where
$\mathbf{u} = (u,v)$
is the velocity,
$h(x,y)$
is the height and motion is on an
$x,y$
plane perpendicular to the vertical unit vector
$\hat {\boldsymbol{z}}$
.
$\mathcal {F}(h) = 1$
for standard shallow water, and
\begin{equation} \mathcal {F}(h) = \frac {1}{h^3} ,\end{equation}
for MSW. The non-dimensional parameters are the Froude and Rossby numbers,
$Fr$
and
$Ro$
, where
$Fr$
is the non-dimensional inverse phase speed of linear gravity waves unaffected by rotation, and
$Ro$
represents the ratio of inertial to Coriolis forces. Throughout, we take
$Fr = 0.3$
and
$Ro = 0.4$
.
The flow is initialised with the output of an incompressible two-dimensional Navier–Stokes simulation in a fully developed turbulent state, with height
$h$
set to be initially in geostrophic balance (as in KV23)
\begin{equation} \frac {1}{Ro}\hat {\boldsymbol {z}}\times \boldsymbol {u} = -\frac {1}{Fr^2}\nabla h, \end{equation}
and is allowed to evolve freely. The non-dimensionalisation of the height is such that
$h = 1 + \eta$
, where the height perturbation
$\eta = 0$
for a flow at rest. For
$\eta \ll 1$
, the MSW term in (4.3) therefore scales as
$\mathcal {F}(h) = 1 + O(\eta )$
, and MSW approximates standard shallow water. For a flow in geostrophic balance,
$\eta \sim Fr^2/Ro$
(see (4.4)), so
$Fr^2/Ro \ll 1$
is the condition for such a geostrophic flow to behave similarly to a standard shallow water flow. Here,
$Fr^2/Ro = 0.225$
, and we find that this is sufficient to prevent any spurious behaviour from the shallow water modification.
We superimpose a linear Poincaré wave on this initial condition and allow it to evolve alongside and interact with the geostrophic turbulence. The linearisation of the MSW (4.1)–(4.2) is identical to the original shallow water system. Linear wave solutions have frequency
$\omega$
satisfying the dispersion relation
\begin{equation} \omega ^2 = \frac {1}{Ro^2} + \frac {|\boldsymbol {k}|^2}{Fr^2}, \end{equation}
where
$\boldsymbol {k} = (k,l)$
is the wavenumber. The height perturbation
$\eta$
of the waves scales as
$A Ro$
, where
$A$
is the maximum amplitude of the vorticity of the initialised wave. We take
$A = 0.5$
, so that
$A Ro = 0.2$
, and the waves are also sufficiently linear to not be obviously affected by the MSW term (4.3), aside from their lack of nonlinear steepening as intended. The mode-1 (
$|\boldsymbol {k}| = 1$
) wave has frequency
$\omega = 4.17$
.
Starting from this initial condition, we evolve the MSW equations alongside the Lagrangian mean equations using a pseudo-spectral solver as in KV23, with a fourth-order Runge–Kutta scheme for the advective terms (Baker et al. Reference Baker, Kafiabad, Maitland-Davies and Vanneste2024). We use a non-dimensional domain size of
$2\unicode{x03C0} \times 2\unicode{x03C0}$
, with 256 gridpoints in the
$x$
and
$y$
directions. A hyperviscous (Laplacian to the power four) term is added to the momentum equation (4.1) to remove energy at small scales. A hyperviscous term can also be added to the Lagrangian mean equations, which is found to be necessary for numerical stability when integrating strategy 1 over long time intervals. However, this is not necessary and is found to introduce error in the resulting Lagrangian mean for strategies 2 and 3, where the forcing terms on the right-hand side of the scalar equations (3.23) and (3.27) seem to stabilise the simulations, so the Lagrangian mean equations of strategies 2 and 3 are run without any viscous terms.
We implement each of strategies 1, 2, and 3, but show results from only strategies 1 and 3, since the results of strategies 2 and 3 are visually identical (although their difference is quantified in § 5), but strategy 3 is faster (see Appendix E). Unless otherwise stated, the equations are run using strategy 3 for an averaging period of
$2T = 40$
, over which time the mean and wave components of the flow both evolve. The weight function
$G$
is truncated to this finite interval, and renormalised to ensure that (2.11) holds exactly over the interval. In the case of a low-pass with weight function given by (2.6)
\begin{equation} \int _{t^{\ast}-T}^{t^{\ast} + T}G_{\mathrm {LP}}(t^{\ast}-s)\,\mathrm {d} s \rightarrow 1 \hspace {0.5cm} \mathrm {as}\ \omega _c T \rightarrow \infty , \end{equation}
so when
$\omega _cT$
is sufficiently large, the normalisation requirement still approximately holds. Unless otherwise stated, we use a low-pass cutoff frequency of
$\omega _c = 2$
, so that
$\omega _c T = 40$
, and
$\int _{t^{\ast}-T}^{t^{\ast} + T}G_{\mathrm {LP}}(t^{\ast}-s)\,\mathrm {d} s = 1.01$
. Appendix F shows the impact of changing
$T$
.
The scalar field to be averaged is the relative vorticity
\begin{equation} \zeta = \frac {\partial v}{\partial x} - \frac {\partial u}{\partial y}. \end{equation}
First, we show results for the Lagrangian mean of the vorticity and compare the different strategies. We then explain the various ways that the flow can be decomposed into wave and mean components, before showing results for these decompositions.
5. Results: Lagrangian mean
Figure 2(a) shows the instantaneous vorticity at the midpoint of the averaging interval for comparison to various means. Whilst there is a high-amplitude mode-1 wave present in the instantaneous vorticity, this wave is removed by the averaging procedures in figures 2(b), 2(c), 2(e) and 2(f).
Figure 2. Shallow water relative vorticity for a simulation over 40 time units (
$T=20$
). The mode-1 wave frequency is
$\omega = 4.17$
, and the low-pass filters use a cutoff frequency of
$\omega _c = 2$
. (a) Instantaneous vorticity at the interval midpoint
$t^{\ast} = 20$
. (b) Lagrangian and (c) Eulerian low-pass at
$t^{\ast} = 20$
. (e) Lagrangian and (f) Eulerian top-hat mean at
$t^{\ast} = 20$
, computed over the interval
$[18,22]$
, i.e.
$T = 2$
. (d) Value of
$G(t)$
for the low-pass and top-hat means, showing that
$T = 2$
is an appropriate averaging interval for the top-hat to compare it to the low-pass. The directory including the Jupyter notebook that generated this figure can be accessed at https://cambridge.org/S0022112025000424/JFM-Notebooks/files/Figure-2/Figure-2.ipynb.
Figures 2(b) and 2(c) show the Lagrangian and Eulerian low-pass means of vorticity. The Lagrangian low-pass retains more of the intensity of the vortices than the Eulerian low-pass, since the effect of the large-amplitude wave displacement on the background turbulence leads to a blurring of the vortices in the Eulerian low-pass. This low-pass is calculated over an interval of 40 time units (
$T = 20$
), and the corresponding weight function
$G(t)$
is shown in figure 2(d). The root-mean-squared difference between the Lagrangian mean vorticity in figure 2(b) when calculated with strategies 2 and 3 is 0.003, with a maximum difference of 0.03.
Also shown in figure 2(d) is
$G(t)$
for a top-hat mean with a comparable averaging time scale of
$T=2$
. Figures 2(e) and 2(f) show the corresponding Lagrangian and Eulerian top-hat means at the same value of
$t^{\ast}$
as for the low-pass. Whilst there are not qualitative differences between the top-hat and low-pass mean vorticity, there are differences that are evident when the fields are viewed as a time series.
Figure 3. As in figure 2, showing the time (
$t^{\ast}$
) evolution of each field at
$y = 2.8$
. The directory including the Jupyter notebook that generated this figure can be accessed at https://cambridge.org/S0022112025000424/JFM-Notebooks/files/Figure-3/Figure-3.ipynb.
Figure 3 shows the same fields as figure 2 over a time series in
$t^{\ast}$
. We note here that figure 3 is not showing the evolution in
$t$
of the partial Lagrangian mean fields described in § 3. Instead, for each value of
$t^{\ast}$
, a set of Lagrangian mean equations are solved over the interval
$[t^{\ast}-T,t^{\ast}+T]$
, and the values of the full Lagrangian means (referenced to time
$t^{\ast}$
) are shown.
Whilst the oscillations at the wave frequency are removed by the Lagrangian low-pass in figure 3(b), they are still evident in the top-hat mean shown in figure 3(e). This is because the top-hat mean is less selective in the frequencies that are filtered. In this simple test case, we could have chosen the averaging period of the top-hat to exactly be the period of the wave to better remove this wave signal. However, in the general case of a continuous spectrum of waves, the top-hat mean would not be able to perfectly remove all waves with a given cutoff frequency.
The Eulerian low-pass mean in figure 3(c) is more effective in removing the wave oscillations than the Lagrangian top-hat, but, as in figure 2(c), the resulting mean flow is blurred. The Eulerian top-hat in figure 3 suffers from both blurring and residual wave signal. Hereafter, we focus on the low-pass filter, as it gives more control over the frequencies to remove, and consider the relative merits of strategies 1 and 3.
Figure 4. Comparison of calculation of
$\overline {\zeta }^{L}$
using strategies 1 and 3 with
$T=20$
: (a)
$\tilde {\zeta }$
, found using strategy 1, (b),
$\overline {\zeta }^{L}$
, found by remapping
$\tilde {\zeta }$
using
$\boldsymbol {\Xi }^{1 \mapsto 2}$
(c), the
$y$
component of
$\boldsymbol {\Xi }^{1 \mapsto 2}$
(d),
$\zeta ^{\ast}$
, found using strategy 3, (e),
$\overline {\zeta }^{L}$
, found by remapping
$\zeta ^{\ast}$
using
$\boldsymbol {\Xi }^{3 \mapsto 2}$
and (f) the
$y$
component of
$\boldsymbol {\Xi }^{3 \mapsto 2}$
. The
$x$
and
$y$
axes correspond to
$x$
and
$y$
coordinates of the full domain. The directory including the Jupyter notebook that generated this figure can be accessed at https://cambridge.org/S0022112025000424/JFM-Notebooks/files/Figure-4/Figure-4.ipynb.
Figure 4(a) shows the direct output
$\tilde {\zeta }$
of strategy 1, and figure 4(d) the direct output
$\zeta ^{\ast}$
of strategy 3 (which are rearrangements of each other). To find
$\overline {\zeta }^{L}$
, these fields are remapped to the trajectory mean coordinate, using
$\boldsymbol {\Xi }^{1\mapsto 2}$
for strategy 1 and
$\boldsymbol {\Xi }^{3\mapsto 2}$
for strategy 3. The
$y$
components of these maps are shown in figures 4(c) and 4(f) respectively, and the resulting
$\overline {\zeta }^{L}$
in figures 4(b) and 4(e). The map
$\boldsymbol {\Xi }^{1\mapsto 2}$
is complicated as it represents the motion of the flow between the trajectory mean and end positions, and the resulting
$\overline {\zeta }^{L}$
is poorly interpolated. However,
$\boldsymbol {\Xi }^{3\mapsto 2}$
differs from the identity only by the wave perturbation, and therefore results in a clean interpolation to
$\overline {\zeta }^{L}$
. Therefore, strategy 3 is preferred over strategy 1 when the mean flow varies significantly over the averaging interval, although in a flow with a more complex ‘wave’ component, the final mapping of strategy 3 may still be too complicated. However, over a shorter averaging interval, where the interpolating map shown in figure 4(c) is simpler, strategy 1 may be more accurate than strategy 3, since strategy 3 can accumulate interpolation errors at each time step. We discuss this further in § 8. The complexity of the strategy 1 mapping in figure 4(c) is also partly due to the long averaging interval used – the use of a filter that is more localised in time than the strict low-pass (such as a Butterworth or Gaussian filter) would allow a shorter averaging interval, and correspondingly a less complex final interpolation in strategy 1 (and also in the forcing terms of the strategy 2 and 3 equations).
Having found the Lagrangian mean of a flow, we now consider how to define the wave-like component of the flow.
6. Lagrangian wave formulation
There are several ways to define the wave-like component of the flow. The first (and perhaps most common) is to define waves as high-frequency perturbations at a fixed point. If the Eulerian mean of some scalar
$f$
is given by
$\overline {f}^{E}(\boldsymbol {x},t)$
, then the Eulerian wave perturbation is defined as
\begin{equation} \overline {f}^{w}_{E}(\boldsymbol {x},t) = f(\boldsymbol {x},t) - \overline {f}^{E}(\boldsymbol {x},t), \end{equation}
where superscript
$\mathrm {w}$
represents the wave component. However, the Lagrangian mean is more effective than the Eulerian mean for recovering a mean flow in the presence of large-amplitude waves (e.g. figures 2 and 3), or when the waves are significantly Doppler shifted (Shakespeare et al. Reference Shakespeare, Gibson, Hogg, Bachman, Keating and Velzeboer2021). We should therefore define waves to be high-frequency motions in the Lagrangian frame, and mean flows to be low-frequency motions in the Lagrangian frame, with some appropriate cutoff frequency separating the two. By this definition, mean flows must not necessarily be balanced in the sense of geostrophic balance, or even slowly varying in the Eulerian reference frame.
However, this definition of a wave is still not precise enough. We have the option to define waves in either a ‘semi-Eulerian’ or a Lagrangian way. We may say a wave perturbation is:
-
(i) Eulerian: an Eulerian high-frequency perturbation at a fixed point.
-
(ii) Semi-Eulerian: a Lagrangian high-frequency perturbation at a fixed point.
-
(iii) Lagrangian: a Lagrangian high-frequency perturbation following a particle.
The difference between these viewpoints stems from the fact that waves impact the flow in two ways: through changing the value of a scalar as seen by a flow-following particle, and through displacing the mean flow. The semi-Eulerian definition of the wave field encompasses these two effects, whereas the Lagrangian wave field only represents the changes in value of a scalar on a particle due to the wave.
Before posing these decompositions mathematically, we consider a simple example to elucidate the difference between them. Consider a two-dimensional (2-D)
$(x,z)$
background flow with uniform velocity
$\textbf {U} = (U_0,0)$
and buoyancy
$B(z)$
that is stably stratified. Steady internal lee waves can propagate on this base state, giving total buoyancy
$b$
and horizontal
$x$
-velocity
$u(x,z)$
of the form
\begin{align} b(x,z) &= B(z) + b^{\prime}(x,z), \end{align}
\begin{align} u(x,z) &= U_0 + u^{\prime}(x,z). \end{align}
Lee waves are generated by flow over topography in the ocean and atmosphere, and are phase locked to topography such that they are steady in the rest frame, hence
$u$
and
$b$
are independent of time. The variables
$b(x,z)$
and
$u(x,z)$
are therefore unchanged by an Eulerian mean, and the Eulerian buoyancy and velocity wave perturbations are zero.
The Lagrangian means of
$b(x,z)$
and
$u(x,z)$
, when filtered with an appropriate cutoff frequency that is lower than the intrinsic wave frequency, are
$B(z)$
and
$U_0$
respectively, thus the semi-Eulerian wave perturbations are
$b^{\prime}$
and
$u^{\prime}$
.
In the absence of diffusion, buoyancy is a conservative tracer satisfying
\begin{equation} \frac {D b}{Dt} = \frac {\partial b}{\partial t} + \boldsymbol {u} \boldsymbol {\cdot } \nabla b = 0, \end{equation}
thus buoyancy is constant following a particle, and the Lagrangian buoyancy wave perturbation is zero. The wave velocity
$u$
is not constant following a particle, therefore the Lagrangian velocity perturbation is non-zero (and unknown for now).
6.1. Semi-Eulerian wave definition
The semi-Eulerian wave field is defined to be the instantaneous field minus the Lagrangian mean field at a fixed spatial location. The mean field is given by the Lagrangian weighted mean
$\overline {f}^{L}$
, which, as discussed in § 2.3, contains no wave signal when weighted with the appropriate frequency filter. We define
\begin{equation} f_{S\mbox{-}E}^{w}(\boldsymbol {x},t^{\ast}) = f(\boldsymbol {x},t^{\ast}) - \overline {f}^{L}(\boldsymbol {x},t^{\ast}), \end{equation}
where the subscript
$S\mbox{-}E$
denotes a semi-Eulerian wave definition.
6.2. Lagrangian wave definitions
The Lagrangian wave field is defined as the Lagrangian high-frequency perturbation on a trajectory. We have two further options for how this is itself defined – either at the Lagrangian trajectory midpoint, or at the Lagrangian mean position
\begin{align} f_{L1}^{w}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast}) &= f(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast}) - f^{\ast}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast}) \end{align}
\begin{align} &= f(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast}) - \overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast})\, , \end{align}
\begin{align} f_{L2}^{w}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast}) &=f(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast}) - \overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast}). \end{align}
Using the map
$\boldsymbol {\Xi }^{3\mapsto 2}$
defined by
$\boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast}) = \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})$
(cf. (3.3)–(3.5)), and letting a dummy variable
$\boldsymbol {x} = \boldsymbol {\varphi }(\boldsymbol {a},t^{\ast})$
in (6.7) and
$\boldsymbol {x} = \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})$
in (6.8) we obtain the alternative forms
\begin{align} f_{L1}^{w}(\boldsymbol {x},t^{\ast}) &= f(\boldsymbol {x},t^{\ast}) - f^{\ast}(\boldsymbol {x},t^{\ast}), \end{align}
\begin{align} &= f(\boldsymbol {x},t^{\ast}) - \overline {f}^{L}(\boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {x},t^{\ast}),t^{\ast})\, , \end{align}
\begin{align} f_{L2}^{w}(\boldsymbol {x},t^{\ast}) &= f((\boldsymbol {\Xi }^{3\mapsto 2})^{-1}(\boldsymbol {x},t^{\ast}),t^{\ast}) - \overline {f}^{L}(\boldsymbol {x},t^{\ast}) .\end{align}
Note that the two definitions (6.10) and (6.11) are just rearrangements of each other such that
$f_{L1}^{w}(\boldsymbol {x},t^{\ast}) = f_{L2}^{w}(\boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {x},t^{\ast}),t^{\ast})$
.
6.3. Comparing wave definitions
Having defined four different wave fields (one Eulerian (6.1), one semi-Eulerian (6.5) and two Lagrangian (6.10) and (6.11)), we now consider the features of each, with a focus on the semi-Eulerian/Lagrangian definitions, having already motivated the Lagrangian over the Eulerian mean. We note that the wave definitions given here do not depend on the strategy with which they are calculated.
The semi-Eulerian definition (6.5) is perhaps the most straightforward to understand, since the wave is defined as ‘what is left when you remove the Lagrangian mean field’. When it is desirable to write the total field as a sum of mean and wave components at the same spatial location, as is done in the lee wave example (6.2) and (6.3), this is the most helpful decomposition. However, although the two terms on the right-hand side of (6.5) are defined at the same spatial location, the instantaneous field
$f$
is evaluated at the position
$\boldsymbol {x}$
which is not necessarily on the path of the particle whose mean position is
$\boldsymbol {x}$
, and whose mean is evaluated in the second term. Hence, we are subtracting the mean of particle from the instantaneous value of a different particle. As in § 2.3, the Lagrangian low-pass filter applied to the wave field would ideally return zero, and this is not the case for the semi-Eulerian wave field when filtered along trajectories of either the original flow
$\boldsymbol {\varphi }$
or the mean flow
$\bar {\bar {\boldsymbol {\varphi }}}$
.
However, the Lagrangian definitions do have this property, and it can be shown that (assuming a simple band-pass filter as described in § 2.3) the first Lagrangian wave field is zero when low-pass filtered along the original flow paths, and the second is zero when low-pass filtered along the mean flow paths. For comparison with the well-known notation of Andrews & McIntyre (Reference Andrews and McIntyre1978), we note that the second Lagrangian wave definition (6.11) corresponds to their Lagrangian disturbance quantities with a superscript
$l$
(e.g. their equation (2.11)), although their Eulerian average (such that
$\overline {\,f^l\,}^{\, E} = 0$
) is replaced in our case with a Lagrangian average along mean flow trajectories since we do not assume separation of time-scales (see Appendix B).
In the first Lagrangian wave definition (6.10), a deformation of the mean field that includes the impact of the wave disturbance of the mean field (
$f^{\ast}$
) is subtracted from the instantaneous field to give a wave component
$f_{L1}^{w}$
that documents only the changes to the value of the field seen by a particle. This is the wave field that is found by filtering methods that use particle tracking with particles seeded at the reference time
$t^{\ast}$
, such as Shakespeare et al. (Reference Shakespeare, Gibson, Hogg, Bachman, Keating and Velzeboer2021), which directly find
$f^{\ast}$
as the mean field (although such methods usually track with horizontal velocities only, so only approximate
$f^{\ast}$
). We will later see that this wave decomposition can give a much clearer view of the wave field than the semi-Eulerian definition, since the wave component does not include the wave-displaced mean flow. However,
$f^{\ast}(\boldsymbol {x},t)$
is not the true mean field, as it includes a wave signal (see figure 6(c)).
The second Lagrangian wave definition (6.11) is similar to the first, but to find the wave field the instantaneous total field must first be deformed to remove the effect of wave displacement, before the mean field is subtracted. Therefore, neither of the Lagrangian descriptions give a decomposition that can be written as wave + mean = total at a fixed spatial location.
Figure 5. The four different wave decompositions: (top) Eulerian, (second row) semi-Eulerian, (third row) Lagrangian first definition and (bottom) Lagrangian second definition. For each row, the middle ‘mean’ field is subtracted from the left ‘instantaneous’ field to give the right ‘wave’ field. The flow parameters are as for figure 2, and strategy 3 is used. The directory including the Jupyter notebook that generated this figure can be accessed at https://cambridge.org/S0022112025000424/JFM-Notebooks/files/Figure-5/Figure-5.ipynb.
7. Results: Lagrangian waves
Figure 5 shows the four different wave decompositions discussed above for the same example as in figure 2, where in each case the left column minus the middle column gives the wave perturbation in the right column. Both the Eulerian and semi-Lagrangian wave definitions give a wave that has a significant signature of the turbulent mean flow. In the Eulerian case, this is because the mean flow is blurred by the high-amplitude wave perturbations, and in the semi-Eulerian case this is because the deformation of the mean field by the high-amplitude wave is included in the wave definition.
The Lagrangian wave definitions in figures 5(i) and 5(l) are much cleaner as they only represent the wave vorticity. However, the original mode-1 plane wave is not perfectly recovered, as can be seen in figures 5(i) and 5(l). This is a result of nonlinear wave–mean interactions. There appear to be two such types of interaction – large-scale deformations of the plane wave due to interaction with the mean flow (seen more clearly in figure 7 and supplementary movie 1), and high-frequency oscillations of the mean flow that appear in figures 5(i) and 5(l) at the same spatial scales as the mean flow. The amplitude of this turbulence-like pattern scales with the wave linearity, is independent of grid resolution (making it unlikely to be caused by interpolation errors – see § 8), and is the same in all strategies. The time evolution of the mean and wave fields is shown in supplementary movie 1.
Figure 6. Hovmöller (space–time) diagrams of vorticity: (a) instantaneous, (b) Lagrangian low-pass, (c) Lagrangian low-pass at the trajectory midpoint, (d) Eulerian low-pass, (e) Lagrangian L1 wave, (f) Lagrangian L2 wave, (g) semi-Eulerian wave and (h) Eulerian wave. Strategy 3 is used to solve for the Lagrangian means at a temporal resolution of 0.2. Parameters are identical to figure 2. All panels are shown at
$y = 2.8$
. The directory including the Jupyter notebook that generated this figure can be accessed at https://cambridge.org/S0022112025000424/JFM-Notebooks/files/Figure-6/Figure-6.ipynb.
Figure 6 shows Hovmöller diagrams of several of the fields shown in figure 5. Comparing the instantaneous (figure 6 a) and Lagrangian mean (figure 6 b) vorticity shows that the wave has been very effectively removed by the Lagrangian low-pass filter.
The wave oscillations are very clear in
$\zeta ^{\ast}$
figure 6(c) – from which
$\overline {\zeta }^{L}$
in figure 6(b) is remapped. Wave oscillations are also visible in the Eulerian mean
$\overline {\zeta }^{\, E}$
in figure 6(d), and the vorticity gradients in the turbulent flow are overly smoothed, as shown in figure 2(c).
The various wave decompositions are shown in the bottom row of figure 6, again demonstrating that the two Lagrangian definitions give a clean representation of the wave field, whereas the Eulerian and semi-Eulerian wave definitions contain significant imprints of the turbulent flow. A movie showing the evolution of
$\zeta, \zeta^{\ast}, \overline{\zeta}^{L}$
and
$\zeta_{L2}^{w}$
over the time series shown in figure 6 is provided in the supplementary material (supplementary movie 1).
Despite the Lagrangian wave perturbations being visually ‘cleaner’ in that they recreate more closely the plane wave with which the simulation was initialised, the physically appropriate wave definition for a given problem is likely to be context-dependent.
Figure 7. An example of different frequency filters with corresponding functions
$\hat {G}(\omega )$
shown in panel (b). (a) Instantaneous vorticity for a MSW simulation with a mode-1 wave in
$x$
and a mode-2 wave in
$y$
of the same vorticity amplitude (
$A = 0.5$
), with respective frequencies 4.17 and 7.12. (c) Lagrangian low-pass filter of the flow in panel a with a cutoff frequency of 2, so that both waves are removed, (d) as in panel (c) with a cutoff frequency of 5.5, so that only the mode-2 wave is removed and (e) as in panel (c) with the filter defined in (2.12) and labelled ‘band-pass’ in panel (b), with
$\omega _1 = 2$
and
$\omega _2 = 5.5$
, so that the mode-2 wave is retained and the mode-1 removed. Panels (f), (g) and (h) show L2 wave perturbation corresponding to panels (c), (d) and (e), respectively. The directory including the Jupyter notebook that generated this figure can be accessed at https://cambridge.org/S0022112025000424/JFM-Notebooks/files/Figure-7/Figure-7.ipynb.
Finally, we present an example to demonstrate the flexibility of Lagrangian frequency filtering. Figure 7(a) shows the instantaneous vorticity of a flow that has been initialised with a turbulent flow and mode-1 wave in the
$x$
direction (as before), and also a mode-2 wave in the
$y$
direction. The waves have the same amplitude in vorticity, and have frequencies 4.17 and 7.12, respectively. In figures 7(c) and 7(f), the mean and L2 wave perturbations are shown for a low-pass filter at cutoff frequency
$\omega _c = 2$
, which removes both waves from the mean flow. Figures 7(d) and 7(g) are as for figures 7(c) and 7(f) with cutoff frequency
$\omega _c = 5.5$
, which retains the mode-1 wave in the ‘mean’ flow (figure 7
d), and leaves the mode-2 wave in the ‘wave’ perturbation (figure 7
g). In figures 7(e) and 7(h) the filter defined in (2.12) is used, such that the Lagrangian mean operation removes frequencies between
$\omega _1 =2$
and
$\omega _2 = 5.5$
and retains all other frequencies. Therefore, the mode-2 wave is kept as part of the ‘mean’, and the mode-1 is in the ‘wave’ perturbation. Each of the weight functions are shown in figure 7(b). As in figure 5, large-scale departures of the wave perturbations in figures 7(g) and 7(h) from a perfect mode-2 and mode-1 plane wave respectively are again attributed to nonlinear interactions between the waves and between the waves and the turbulent flow.
8. Numerical errors and interpolation
There are two primary error sources in our method – the truncation of the interval length and the interpolations, and these error sources are shared by particle tracking methods. Appendix F shows the impact of varying the half-interval length
$T$
. We found that increasing
$T$
reduced the remaining wave-frequency oscillations, but increasing the interval length past
$T = 20$
made negligible difference to the solutions, making
$T=20$
our value of choice (so that
$\omega _c T = 40$
, and the averaging interval
$2T$
is 12.7 times longer than the cutoff period). It may, however, be worth reducing
$T$
and suffering a small error of this type to reduce the computational expense of the calculations.
The other source of error comes from interpolation. In the methods discussed here, there are two types of interpolation. The first is performed at every time step of strategies 2 and 3 on the right-hand sides of (3.20), (3.23), (3.27) and (3.28), and requires finding a scalar field (the scalar to be averaged, or each component of velocity) at some coordinate
$\boldsymbol {\Xi }^{2\mapsto 1}_p(x,t)$
(strategy 2) or
$\boldsymbol {\Xi }^{3\mapsto 1}_p(x,t)$
(strategy 3), where the scalar is known on a regular grid. The second is the final remapping of
$\tilde {f}$
(in strategy 1) or
$f^{\ast}$
(in strategy 3) to
$\overline {f}^{L}$
. In this case,
$\overline {f}^{L}$
is known at the irregularly spaced locations
$\boldsymbol {\Xi }^{1\mapsto 2}(\boldsymbol {x},t^{\ast})$
or
$\boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {x},t^{\ast})$
, and needs to be found at regularly gridded locations. We do not expect this second type of interpolation to be problematic in strategy 3, as demonstrated in figure 4.
When the flow is such that points advected by the flow become very far apart or very close together over the interval of interest, as represented by an interpolating map with sharp gradients (e.g.
$\boldsymbol {\Xi }^{1\mapsto 2}$
shown in figure 4
c), both types of interpolation can be prone to error. This can be the case when the flow is compressible (or equivalently when performed on a 2-D surface in a 3-D incompressible flow; Shakespeare et al. Reference Shakespeare, Gibson, Hogg, Bachman, Keating and Velzeboer2021), but also when the flow is incompressible and straining or shearing.
In the shallow water case here, the flow is both compressible and straining/shearing, and the mean flow evolves significantly over the time interval (e.g. figure 3 b). We therefore expect to accumulate the first type of interpolation error at each time step, although we do not see evidence of this error in our experiments. The mean flow ‘imprint’ in the wave component in figures 5(i), 5(l), 6(c) and 6(f) is independent of resolution and scales with wave nonlinearity, so we attribute this to physical nonlinear interactions as discussed in § 7.
Although strategies 2 and 3 accumulate interpolation error at every time step, the interpolation terms in the Lagrangian scalar equations (3.23) and (3.27) are weighted by
$G(t^{\ast} - t)$
, which becomes small as the time moves away from the interval midpoint (see figure 2
d). Therefore, the (more accurate) interpolation along trajectories over short times is more important, and long time interpolations become negligible. Increasing the interval time to reduce truncation error does not significantly increase interpolation error in strategy 3.
Strategy 1 does not suffer from accumulation of interpolation error as particle tracking and the other strategies do, needing only one interpolation, but this interpolation can be complex and inaccurate (as shown in figure 4), and depends strongly on the length of the averaging interval. Thus increasing the averaging interval worsens interpolation error in strategy 1, but improves truncation error.
A choice over whether strategy 1, 2 or 3 is optimal should be made based on the nature of the flow, time scales, boundary conditions, computational parallelisation and weight function involved.
9. Discussion
In this work, we have extended the PDE-based approach of KV23 for finding a top-hat Lagrangian mean to a Lagrangian mean with a general convolutional weight function. In particular, this has allowed us to present a method for Lagrangian frequency filtering, whereby specific intrinsic frequencies of a flow can be isolated from the rest of the flow. We have also derived some of the special properties of Lagrangian mean flows that hold for particular weight functions, and explored several different wave–mean decompositions.
In addition to re-deriving the strategies 1 and 2 of KV23 for a general weight function, we have presented a novel strategy 3 that removes some of the difficulties associated with strategies 1 and 2, and have shown that this strategy allows a clean decomposition of geostrophic turbulence and large-amplitude Poincaré waves in a simple rotating shallow water system.
In the system presented here, Lagrangian filtering aims to recover the mean flow without the signature of the large-amplitude wave displacements. We have demonstrated the ability of our method to achieve this. However, an equally important use of Lagrangian filtering is to allow decomposition of waves and mean flows when the waves have been Doppler shifted by the mean flow, that is when the flow speed is large compared to the phase speed of the waves. This is the use that Shakespeare et al. (Reference Shakespeare, Gibson, Hogg, Bachman, Keating and Velzeboer2021) focus on in their presentation of Lagrangian filtering. Although we do not show an example of Doppler shifting of waves by the mean flow here, the method is straightforwardly applicable.
Our method for Lagrangian filtering can be compared to existing particle tracking methods. Lagrangian filtered fields can be found by tracking particles online, although Lagrangian means are then defined at and remapped from the initially seeded particle positions. This can lead to problems with particle clustering similar to those discussed in § 8. To tackle this problem, Shakespeare et al. (Reference Shakespeare, Gibson, Hogg, Bachman, Keating and Velzeboer2021) recently presented an open-source implementation of offline Lagrangian filtering. Their method uses offline simulation data (scalars and horizontal velocities) to track particles backwards and forwards from the interval midpoint, finding the time series of a scalar on a particle, temporally filtering and assigning the filtered scalar to the trajectory midpoint. This method (when carried out using 3-D velocities) directly finds
$f^{\ast}$
(in our notation), rather than the generalised Lagrangian mean
$\overline {f}^{L}$
, but this may be sufficient if the waves are low amplitude and
$f^{\ast} \simeq \overline {f}^{L}$
, or if the
$L1$
wave decomposition is needed. Alternatively,
$\overline {f}^{L}$
could be recovered by finding the Lagrangian mean of position and performing an interpolation of
$f^{\ast}$
to
$\overline {f}^{L}$
(similar to the final step in our strategy 3). Offline particle tracking requires saving, storing and processing large quantities of simulation output, requiring high storage and post-processing cost. Particle tracking also suffers from expense and error associated with interpolation at each time step, similarly to our strategies 2 and 3.
In contrast, our method solves the Lagrangian mean equations at the same time as the evolution equations of the flow itself, so saving high-frequency simulation output is not required. This also allows the Lagrangian equations to be solved on the same grid and using the same numerical scheme as the original simulation. There is flexibility over the weight function used and the specific Lagrangian mean and wave definitions that are solved for. However, this does increase the computational expense of the simulation itself. In our 2-D shallow water example, using (the cheapest) strategy 1 to solve for
$\overline {f}^{L}$
increased the computation time over shallow water alone by 66 %, and using strategy 3 by 132 % (see Appendix E).
We therefore expect the two different methods to have different uses. When filtering existing simulation output, or output from a large and complex general circulation model, it may be preferable to use particle tracking offline. However, for process studies where wave or mean identification is a primary objective, our method is easily implemented, more flexible and requires much less storage.
When finding the Lagrangian mean at high temporal resolution, the expense of our method increases greatly since a set of Lagrangian mean equations needs to be solved for each time
$t^{\ast}$
where the Lagrangian mean is required (e.g. as in figure 6). Particle tracking methods also suffer from this drawback to some extent. If a ‘slow’ Lagrangian mean is the quantity of interest, then a time series of Lagrangian means can be found at a coarse time resolution that only resolves this slow variation. The semi-Eulerian wave perturbation can then be found at all times by interpolating this mean to the time of interest and removing it from the instantaneous field. However, if one of the Lagrangian wave definitions is needed, then the Lagrangian mean calculation needs to be carried out at a temporal resolution that captures the waves.
There is, however, a special class of weight functions that are exponential or sum-of-exponentials, in which case the partial Lagrangian mean found at each time step during the evolution of the Lagrangian mean equations is itself the full mean (Minz et al. Reference Minz, Baker, Kafiabad and Vanneste2024). This allows the Lagrangian mean to be found at each time step with the expense of solving only one set of Lagrangian mean equations. The drawback of this method is an inability to freely choose the weight function (for example, to filter at a specific frequency) but the significant improvement in computational cost may make this worthwhile. The derivation and an evaluation of the exponential mean is presented in Minz et al. (Reference Minz, Baker, Kafiabad and Vanneste2024).
We presented three separate strategies for finding the Lagrangian mean, each of which has its own advantages and disadvantages. Strategy 1 is cheapest and is simple to implement (particularly in distributed memory parallelisation), but fails when the mean flow varies significantly over the averaging interval (e.g. figure 4 b). Strategy 2 directly finds the Lagrangian mean, but cannot easily be used in bounded domains and is the slowest method (see Appendix E).
Strategy 3 is slower than strategy 1, but faster than strategy 2. It has simple boundary conditions, and can be solved in a periodic or bounded domain. Like strategy 1, it requires a final interpolation step, but this interpolation is simpler than that in strategy 1, and likely to be accurate for a low-pass filter. Future work will focus on implementing the three strategies in a 3-D Boussinesq solver and testing their ability for Lagrangian filtering of different flow configurations.
Supplementary movie
Supplementary movie and Computational Notebook files are available at https://doi.org/10.1017/jfm.2025.42. Computational Notebooks can also be found online at https://www.cambridge.org/S0022112025000424/JFM-Notebooks.
Acknowledgements
We are grateful to the editor and to two anonymous reviewers for their helpful and constructive comments.
Funding
L.E.B. is supported by the Engineering and Physical Sciences Research Council, grant EP/X028135/1. H.A.K. is supported by the Engineering and Physical Sciences Research Council, grant EP/Y021479/1. J.V. is supported by the UK Natural Environment Research Council, grant NE/W002876/1.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available at https://doi.org/10.5281/zenodo.14237745 (Baker et al. Reference Baker, Kafiabad, Maitland-Davies and Vanneste2024).
Appendix A. The mean of a Lagrangian mean scalar
As explained in § 2.3, we would like the mean
$\overline {f}^{L}(\boldsymbol {x},t^{\ast})$
to satisfy a property that we expect of a mean, namely that the mean is unchanged by reapplying the averaging operation (2.13).
Here, we introduce some clarifying notation to define what is meant by the operation on the left-hand-side of (2.13), since the averaging operation itself depends on the flow with respect to which the Lagrangian mean is taken. We define Lagrangian averaging operators that act on some scalar field
$h(\boldsymbol {x},t)$
\begin{align} \overline {\,h\,}^{\,\boldsymbol {\varphi }}(\boldsymbol {a},t^{\ast}) &= \int _{-\infty }^{\infty } G(t^{\ast}-s)h(\boldsymbol {\varphi }(\boldsymbol {a},s),s) \,\mathrm {d} s ,\end{align}
\begin{align} \overline {\,h\,}^{\,\bar {\bar {\boldsymbol {\varphi }}}}(\boldsymbol {a},t^{\ast}) &= \int _{-\infty }^{\infty } G(t^{\ast}-s)h(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},s),s) \,\mathrm {d} s, \end{align}
such that
$\overline {\,(\cdot )\,}^{\,\boldsymbol {\varphi }}$
denotes a Lagrangian mean at time
$t^{\ast}$
along the flow defined by
$\boldsymbol {\varphi }$
for a trajectory labelled by
$\boldsymbol {a}$
, and
$\overline {\,(\cdot )\,}^{\,\bar {\bar {\boldsymbol {\varphi }}}}$
similar for a (mean) flow defined by
$\bar {\bar {\boldsymbol {\varphi }}}$
. The definition of the generalised Lagrangian mean
$\overline {f}^{L}$
of a scalar
$f$
given in (2.9) can then be written
\begin{equation} \overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast}) \equiv \left (\overline {f}^{L}\circ \bar {\bar {\boldsymbol {\varphi }}}\right ) (\boldsymbol {a},t^{\ast}) \equiv \overline {\,f\,}^{\,\boldsymbol {\varphi }}(\boldsymbol {a},t^{\ast}), \end{equation}
where function composition is denoted by ‘
$\circ$
’ and taken to apply to the first argument of a given function. Note that
$\overline {\,f\,}^{\,\boldsymbol {\varphi }}$
is a function of the label space
$(\boldsymbol {a},t^{\ast})$
, making it a fully Lagrangian variable, whereas
$\overline {f}^{L}$
is a function of physical space (specifically, the mean position).
The Lagrangian mean of the Lagrangian mean scalar should be taken with respect to the mean flow, so (2.13) can now be posed more carefully as
\begin{equation} \overline {\,\overline {f}^{L}\,}^{\,\bar {\bar {\boldsymbol {\varphi }}}} = \overline {\,f\,}^{\,\boldsymbol {\varphi }}. \end{equation}
We find the conditions on the weight function
$G$
for which (A4) holds. We have
\begin{align} \overline {\,\overline {f}^{L}\,}^{\,\bar {\bar {\boldsymbol {\varphi }}}}(\boldsymbol {a},t^{\ast}) &= \int _{-\infty }^{\infty } G(t^{\ast} -s)\overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},s),s)\,\mathrm {d} s \nonumber \\ &= \int _{-\infty }^{\infty }\left [\int _{-\infty }^{\infty } G(t^{\ast} - u)G(u-s)\,\mathrm {d} u\right ] f(\boldsymbol {\varphi }(\boldsymbol {a},s),s)\,\mathrm {d} s. \end{align}
By comparison with (2.9), we see that
\begin{align}\overline {\,\overline {f}^{L}\,}^{\bar {\bar {\boldsymbol {\varphi }}}}(\boldsymbol {a},t^{\ast}) = \overline {\,f\,}^{\,\boldsymbol {\varphi }}(\boldsymbol {a},t^{\ast}) &\Leftrightarrow \int _{-\infty }^{\infty } G(t^{\ast} - u)G(u-s)\,\mathrm {d} u = G(t^{\ast}-s) \nonumber \\ &\Leftrightarrow \left (\hat {G}(\omega )\right )^2 = \hat {G}(\omega ), \end{align}
where
$\hat {G}(\omega )$
is the Fourier transform of
$G$
, defined in (2.5).
From (A6), we see that condition (A4) is only satisfied if
$\hat {G}(\omega ) = 0$
or
$1$
, or some piece-wise combination of each.
Appendix B. Relation to classical GLM theory and time-scale separation
The development of GLM theory by Andrews & McIntyre (Reference Andrews and McIntyre1978) defines averaging procedures in an abstract way that apply similarly to spatial, temporal or ensemble Lagrangian averages. However, for a time average it is assumed that there is a time-scale separation between the ‘slow’ and ‘fast’ motions to be separated by Lagrangian averaging (Bühler Reference Bühler2014). Here, we explain this assumption and how it relates to our formulation.
First, we define the lift map
$\boldsymbol {\Xi }$
from the mean flow map to the flow map (this is equivalent to
$\boldsymbol {\Xi }^{2 \mapsto 3}$
in the notation of the current work) by
\begin{equation} \boldsymbol {\Xi }(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast}) = \boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}), \end{equation}
so that a particle at position
$\boldsymbol {\Xi }(\boldsymbol {x},t)$
has mean position
$\boldsymbol {x}$
at time
$t$
. In our notation (using the definitions (A1)–(A2)), we have
\begin{equation} \overline {\,f \circ \boldsymbol {\Xi }\,}^{\,\bar {\bar {\boldsymbol {\varphi }}}} = \overline {\,f\,}^{\,\boldsymbol {\varphi }}. \end{equation}
Andrews & McIntyre (Reference Andrews and McIntyre1978) define the Lagrangian mean (which we distinguish from our definition by a prime) as
\begin{equation} {\overline {f}}^{L^{\prime}}(\boldsymbol {x},t^{\ast}) = \overline {\,f \circ \boldsymbol {\Xi }\,}^{\, E}(\boldsymbol {x},t^{\ast}), \end{equation}
where
$\overline {(\cdot)}^{E}$
denotes the Eulerian time average defined in (2.10).
Rewriting our definition of the Lagrangian mean for comparison with (B3) using (A3) and (B2) gives
\begin{equation} \overline {f}^{L} \circ \bar {\bar {\boldsymbol {\varphi }}} = \overline {\,f \circ \boldsymbol {\Xi }\,}^{\,\bar {\bar {\boldsymbol {\varphi }}}}. \end{equation}
Therefore, the definitions of
$\overline {f}^{L}$
and
$\overline {f}^{L^{\prime}}$
are equivalent (and the assumption of separate slow and fast time scales holds) when
\begin{equation} \overline {\,f \circ \boldsymbol {\Xi }\,}^{E}\circ \bar {\bar {\boldsymbol {\varphi }}} = \overline {\,f \circ \boldsymbol {\Xi }\,}^{\bar {\bar {\boldsymbol {\varphi }}}}, \end{equation}
or equivalently, when the Eulerian mean
$\overline {\,(\cdot )\,}^{\, E}$
in definition (B3) is a good approximation to the more general expression
$\overline {\,(\cdot )\,}^{\,\bar {\bar {\boldsymbol {\varphi }}}}\circ \bar {\bar {\boldsymbol {\varphi }}}^{-1}$
.
To demonstrate this condition in a different way, we can write
\begin{align} \overline {f}^{L^{\prime}}(\boldsymbol {x},t^{\ast}) &= \overline {\,f \circ \boldsymbol {\Xi }\,}^{\, E}(\boldsymbol {x},t^{\ast}) \nonumber\\ &= \int _{-\infty }^{\infty } G(t^{\ast}-s) f(\boldsymbol {\Xi }(\boldsymbol {x},s),s) \,\mathrm {d} s. \end{align}
Letting
$\boldsymbol {x} = \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})$
gives
\begin{equation} \overline {f}^{L^{\prime}}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast}) = \int _{-\infty }^{\infty } G(t^{\ast}-s)f(\boldsymbol {\Xi }(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),s),s)\,\mathrm {d} s. \end{equation}
However, by our definitions (A1) and (A2),
\begin{align} \overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast})&= \int _{-\infty }^{\infty } G(t^{\ast}-s)f(\boldsymbol {\varphi }(\boldsymbol {a},s),s),s)\,\mathrm {d} s \nonumber \\ &= \int _{-\infty }^{\infty } G(t^{\ast}-s)f(\boldsymbol {\Xi }(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},s),s),s)\,\mathrm {d} s \, . \end{align}
The two definitions
$\overline {f}^{L}$
and
$\overline {f}^{L^{\prime}}$
are approximately equal provided that (comparing (B7) and (B8))
$f(\boldsymbol {\Xi }(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),s),s) \approx f(\boldsymbol {\Xi }(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},s),s),s)$
where
$G(t^{\ast}-s)$
is not small. This assumes that the mean flow is ‘frozen’ during the averaging operation and is the implicit assumption of time-scale separation behind the Andrews & McIntyre (Reference Andrews and McIntyre1978) definition of the Lagrangian mean.
A flow which illustrates the difference between the two formulations is the classical lee-wave problem discussed in § 6. Consider the flow defined in (6.2) and (6.3). The flow is steady, so
$\boldsymbol {\Xi }$
(as defined in (B1)) and the scalar
$f$
to be averaged are independent of time. Then we have
\begin{align} \overline {f}^{L^{\prime}}(\boldsymbol {x}) &= \int _{-\infty }^{\infty } G(t^{\ast}-s) f(\boldsymbol {\Xi }(\boldsymbol {x})) \,\mathrm {d} s \nonumber \\ &= f(\boldsymbol {\Xi }(\boldsymbol {x})), \end{align}
whereas
\begin{equation} \overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast})) = \int _{-\infty }^{\infty } G(t^{\ast}-s)f(\boldsymbol {\varphi }(\boldsymbol {a},s))\,\mathrm {d} s. \end{equation}
Suppose the scalar
$f$
is the vertical velocity
$w$
. The first definition (B9) is simply a rearrangement of the scalar field and no averaging is performed, so the Lagrangian mean will be non-zero. However, the second definition (as used in this paper) will average over the oscillations of
$w$
on a particle and give zero Lagrangian mean (when
$G$
is defined to remove wave frequencies) as is expected. The difference between the two definitions results from relaxation of the assumption in the second that the mean flow is ‘frozen’ during the averaging interval.
Appendix C. The Lagrangian material derivative
We show that, for convolutional weight functions
$G(t^{\ast}-s)$
, there is a powerful relation between the material derivative of a Lagrangian mean quantity and the Lagrangian mean of a material derivative, namely (cf. (2.16) and (2.17), repeated here)
\begin{equation} \bar {\bar {D}}\overline {f}^{L} = \overline {\,D f\,}^{L}, \end{equation}
where
\begin{align} D &\equiv \frac {\partial }{\partial t} + \boldsymbol {u}\boldsymbol {\cdot } \nabla , \end{align}
\begin{align} \bar {\bar {D}} &\equiv \frac {\partial }{\partial t^{\ast}} + \bar {\bar {\boldsymbol {u}}}\boldsymbol {\cdot } \nabla . \end{align}
We have
\begin{align} \left (\bar {\bar {D}}\overline {f}^{L}\right ) \circ \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}) &= \left ( \frac {\partial \overline {f}^{L}}{\partial t^{\ast}} + \bar {\bar {\boldsymbol {u}}}\boldsymbol {\cdot } \nabla \overline {f}^{L} \right ) \circ \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}) \end{align}
\begin{align} &= \frac {\mathrm {d}}{\mathrm {d} t^{\ast}}\left (\overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast})\right ) \end{align}
\begin{align} &= \frac {\mathrm {d}}{\mathrm {d} t^{\ast}}\int _{-\infty }^\infty G(t^{\ast}-s)f(\boldsymbol {\varphi }(\boldsymbol {a},s),s)\,\mathrm {d} s \end{align}
\begin{align} &= \int _{-\infty }^\infty G^{\prime}(t^{\ast}-s) f(\boldsymbol {\varphi }(\boldsymbol {a},s),s)\,\mathrm {d} s\, \end{align}
\begin{align} &= \int _{-\infty }^\infty G(t^{\ast}-s) \frac {\mathrm {d}}{\mathrm {d} s}\left (f(\boldsymbol {\varphi }(\boldsymbol {a},s),s)\right ) \, \mathrm {d} s \end{align}
\begin{align} &= \int _{-\infty }^\infty G(t^{\ast}-s)\left (\frac {\partial f}{\partial s}(\boldsymbol {\varphi }(\boldsymbol {a},s),s) + \frac {\partial \boldsymbol {\varphi }}{\partial s}(\boldsymbol {a},s) \boldsymbol {\cdot } \nabla f(\boldsymbol {\varphi }(\boldsymbol {a},s),s)\right )\, \mathrm {d} s \end{align}
\begin{align} & = \left (\overline {\,D f\,}^{L}\right )\circ \bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}), \end{align}
where from (C4) and (C5) we used the definition (2.15) of Lagrangian mean velocity, and from (C7) and (C8) we relied on the convolutional form of the weight function and used integration by parts, assuming from (2.11) that
$G(t^{\ast}-s) \rightarrow 0$
as
$s\rightarrow \pm \infty$
. Equivalently, it can be shown that (C1) only holds when the weight function takes the specific (convolutional) form of a frequency filter.
Appendix D. Derivation of strategy 3
Here we derive strategy 3, which solves directly for
$f^{\ast}(\boldsymbol {\varphi }(\boldsymbol {a},t^{\ast}),t^{\ast})$
. For this we need to solve for the map
$\boldsymbol {\Xi }^{3\mapsto 1}$
, and also find a map
$\boldsymbol {\Xi }^{3\mapsto 2}$
to enable us to find
$\overline {f}^{L}(\bar {\bar {\boldsymbol {\varphi }}}(\boldsymbol {a},t^{\ast}),t^{\ast})$
.
We first consider the case
$t \lt t^{\ast}$
. This case is equivalent to strategy 1, since we solve at
${\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t) = \boldsymbol {\varphi }(\boldsymbol {a},t)$
. Differentiating (3.7), or by comparison with (3.12) in strategy 1, we can write an equation for the partial mean
$f^{\ast}_p$
\begin{equation} \frac {\partial f^{\ast}_p}{\partial t}(\boldsymbol {x},t) + \boldsymbol {u}(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla f^{\ast}_p(\boldsymbol {x},t) = f(\boldsymbol {x},t)G(t^{\ast} - t). \end{equation}
We note that for
$t \lt t^{\ast}$
,
$\boldsymbol {\Xi }^{3\mapsto 1}_p(\boldsymbol {x},t) = \boldsymbol {x}$
is the identity map, so, defining
\begin{equation} \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t) = \boldsymbol {\Xi }^{3\mapsto 1}_p(\boldsymbol {x},t) - \boldsymbol {x}, \end{equation}
we have
\begin{equation} \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t) = 0\,, \hspace {1cm} t\lt t^{\ast}. \end{equation}
We also solve for
$\boldsymbol {\Xi }^{3\mapsto 2}(\boldsymbol {x},t^{\ast})$
to enable us to find
$\overline {f}^{L}(\boldsymbol {x},t^{\ast})$
. To do this, we differentiate the definition of
$\boldsymbol {\Xi }^{3\mapsto 2}_p$
in (3.11). For
$t \lt t^{\ast}$
, this is equivalent to solving for
$\boldsymbol {\Xi }^{1\mapsto 2}_p$
in strategy 1, so we define (cf. (3.14))
\begin{equation} \boldsymbol {\xi }^{3\mapsto 2}_p(\boldsymbol {x},t) = \boldsymbol {\Xi }^{3\mapsto 2}_p(\boldsymbol {x},t) - \boldsymbol {x}, \end{equation}
then from (3.13)
\begin{equation} \frac {\partial \boldsymbol {\xi }^{3\mapsto 2}_p}{\partial t}(\boldsymbol {x},t) + \boldsymbol {u}(\boldsymbol {x},t)\boldsymbol {\cdot } \nabla \boldsymbol {\xi }^{3\mapsto 2}_p (\boldsymbol {x},t) = - \boldsymbol {u}(\boldsymbol {x},t)\int _{t^{\ast}-T}^t G(t^{\ast}-s)\,\mathrm {d} s. \end{equation}
We now consider the case
$t \gt t^{\ast}$
. Differentiating (3.7) with respect to
$t$
, we find
\begin{equation} \frac {\partial f^{\ast}_p}{\partial t}(\boldsymbol {x},t) = G(t^{\ast} - t)f(\boldsymbol {x} + \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t)). \end{equation}
We now find an equation for
$\boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t)$
. Differentiating the definition of
$\boldsymbol {\Xi }^{3\mapsto 1}_p$
in (3.10) gives
\begin{align} \frac {\partial \boldsymbol {\Xi }^{3\mapsto 1}_p}{\partial t} &= \boldsymbol {u}(\boldsymbol {\Xi }^{3\mapsto 1}_p(\boldsymbol {x},t),t) \nonumber\\ \Leftrightarrow \hspace {1cm} \frac {\partial \boldsymbol {\xi }^{3\mapsto 1}_p}{\partial t} &= \boldsymbol {u}(\boldsymbol {x} + \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t),t). \end{align}
Finally, we find an equation for
$\boldsymbol {\xi }^{3\mapsto 2}_p(\boldsymbol {x},t)$
. Differentiating the definition of
$\boldsymbol {\Xi }^{3\mapsto 2}_p$
in (3.11) gives
\begin{align} \frac {\partial \boldsymbol {\Xi }^{3\mapsto 2}_p}{\partial t}({\boldsymbol {\varphi }}^{\ast}_p(\boldsymbol {a},t),t) &= \boldsymbol {u}(\boldsymbol {\varphi }(\boldsymbol {a},t),t)\left (1 - \int _{t^{\ast}-T}^t G(t^{\ast}-s)\,\mathrm {d} s\right )\nonumber\\ \Leftrightarrow \hspace {1cm} \frac {\partial \boldsymbol {\xi }^{3\mapsto 2}_p}{\partial t}(\boldsymbol {x},t) &= \boldsymbol {u}(\boldsymbol {x} + \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t),t)\left (1 - \int _{t^{\ast}-T}^t G(t^{\ast}-s)\,\mathrm {d} s\right ). \end{align}
The final equations are summarised as (3.27)–(3.29) in the main text.
Appendix E. Run time of each strategy
Table 2 presents the run time of strategies 1, 2 and 3 when solving for different combinations of the Lagrangian fields that may be required. Simulations are run at
$256 \times 256$
horizontal resolution over 20 time units (
$T=10$
) with a time step of 0.005 (4000 time steps), and an average time taken over three runs. The time reported is for the simulation only – the final remapping in each case takes the same time as
${\lt } 100$
time steps. Vorticity is the only Lagrangian mean scalar being solved for.
Table 2. Run times for the shallow water simulation and Lagrangian mean computation for each strategy, using the code given in Baker et al. (Reference Baker, Kafiabad, Maitland-Davies and Vanneste2024). Times are normalised by the time taken for strategy 1 when solving for
$\overline {f}^{L}$
only. For comparison, when the simulation is run without the Lagrangian mean equations (shallow water only), the corresponding normalised time is 0.6.
The variation in run times between each column of table 2 for each strategy is due to the number of equations being solved and the complexity of these equations. The combinations of PDEs needed for each strategy are listed in table 1.
For each combination of Lagrangian fields to be solved, strategy 1 is the fastest because there is no interpolation on the right-hand side of the Lagrangian mean equations (3.12), (3.13) and (3.16). However, this can come at the expense of accuracy due to the final interpolation needed to recover
$\overline {f}^{L}$
(see figure 4).
The most expensive operations in this pseudo-spectral solver are finding interpolations and calculating nonlinear terms. Strategy 3 is faster than strategy 2 because strategy 2 requires finding both a nonlinear advection term and an interpolated term at every time step in (3.20), (3.23) and (3.24), whereas strategy 3 ((3.27), (3.28) and (3.29)) requires either computing nonlinear terms or computing interpolations at each time step, not both (since
$\boldsymbol {x} + \boldsymbol {\xi }^{3\mapsto 1}_p(\boldsymbol {x},t) = \boldsymbol {x}$
for
$t \lt t^{\ast}$
from (D3)).
Figure 8. Hovmöller diagrams of vorticity for a range of averaging interval times. (a) Instantaneous
$\zeta$
, (b–f) Lagrangian mean
$\overline {\zeta }^{L}$
and (g–k) L2 wave
$\zeta _{L2}^{w}$
. The directory including the Jupyter notebook that generated this figure can be accessed at https://cambridge.org/S0022112025000424/JFM-Notebooks/files/Figure-8/Figure-8.ipynb.
Appendix F. Comparison of filter interval times
Throughout this study, we used an averaging interval time of
$2T = 40$
. Longer averaging times improve the wave decomposition, since there is less truncation error when approximating the full Lagrangian mean (2.9) by the integral over the finite interval (3.2). The condition for the truncation to approximate the full interval is
$\omega _c T \gg 1$
, where
$\omega _c = 2$
here.
Figure 8 shows the impact of increasing the interval time
$2T$
on the time series of Lagrangian mean and L2 wave perturbation. As
$T$
increases, the quality of the filter improves and progressively more wave signal is removed from the Lagrangian mean. The error decreases until
$2T = 40$
. Filters that are more localised in time (such as a Butterworth or Gaussian filter) would also allow earlier truncation and a shorter averaging interval.
Open access
