# Computing Lagrangian means

Hossein A. Kafiabad1<EMAIL_ADDRESS>Jacques Vanneste1  1School of
Mathematics and Maxwell Institute for Mathematical Sciences, University of
Edinburgh, EH9 3FD, UK

###### Abstract

Lagrangian averaging plays an important role in the analysis of wave–mean-flow
interactions and other multiscale fluid phenomena. The numerical computation
of Lagrangian means, e.g. from simulation data, is however challenging.
Typical implementations require tracking a large number of particles to
construct Lagrangian time series which are then averaged using a low-pass
filter. This has drawbacks that include large memory demands, particle
clustering and complications of parallelisation. We develop a novel approach
in which the Lagrangian means of various fields (including particle positions)
are computed by solving partial differential equations (PDEs) that are
integrated over successive averaging time intervals. We propose two
strategies, distinguished by their spatial independent variables. The first,
which generalises the algorithm of Kafiabad (2022, J. Fluid Mech. 940, A2),
uses end-of-interval particle positions; the second directly uses the
Lagrangian mean positions. The PDEs can be discretised in a variety of ways,
e.g. using the same discretisation as that employed for the governing
dynamical equations, and solved on-the-fly to minimise the memory footprint.
We illustrate the new approach with a pseudospectral implementation for the
rotating shallow-water model. Two applications to flows that combine vortical
turbulence and Poincaré waves demonstrate the superiority of Lagrangian
averaging over Eulerian averaging for wave–vortex separation.

##  1 Introduction

Time averaging is a basic yet essential tool in fluid dynamics because of the
ubiquity of phenomena involving multiple time scales. Atmospheric and oceanic
flows, for instance, can often be decomposed usefully into fast and slow
components (e.g. Vallis, 2017; Vanneste, 2013). Time averaging is then used to
interpret numerical-simulation and observational data and to assimilate the
latter in ocean and weather models. More broadly, it can be argued that all
numerical simulations of time-dependent flows rely implicitly on time
averaging, since they do not resolve phenomena on time scales shorter than the
time step.

Time averaging can be performed in different ways. The most straightforward
way is to average the time series of flow variables at fixed spatial positions
to obtain the so-called Eulerian mean. An alternative is to average flow
variables along particle trajectories, that is, at fixed particle labels
instead of fixed positions, to obtain the Lagrangian mean. There are several
practical and conceptual reasons to view Lagrangian averaging as superior to
Eulerian averaging.

Figure 1: Vorticity field and its Lagrangian and Eulerian means in the
shallow-water simulation of §4.2.1: (a) instantaneous vorticity, (b)
Lagrangian mean vorticity, and (c) Eulerian mean vorticity. Both mean fields
share the same averaging period corresponding to 3.6 wave periods. Panels (b)
and (c) share the same colourbar shown to the right of panel (c).

From a practical viewpoint, the advection of fast motion (often consisting of
waves) by slowly evolving flow and vice versa adversely affect their Eulerian
decomposition. For instance, a strong, slowly evolving flow ‘Doppler shifts’
the frequency of fast waves. This can lead to wave frequencies observed at
fixed spatial locations that are much smaller than the intrinsic frequency. It
therefore obscures the time-scale separation between waves and flow.
Lagrangian averaging resolves this issue, as demonstrated in several studies
(Nagai et al., 2015; Shakespeare & Hogg, 2017; Bachman et al., 2020;
Shakespeare et al., 2021; Jones et al., 2022). Conversely, the advection of a
slow background flow by waves leads to flow features that are blurred by
Eulerian averaging but not by Lagrangian averaging. We illustrate this in
figure 1 showing the vorticity field in a simulation of a turbulent rotating
shallow-water flow interacting with a mode-1 Poincaré wave. The instantaneous
vorticity in panel (a) is dominated by the high-amplitude wave. Both the
Lagrangian and Eulerian means (panels (b) and (c) respectively) filter out the
wave, but the Eulerian mean blurs out fine vorticity structures which are well
resolved by the Lagrangian mean (details of this simulation are presented in
§4.2.1).

From a conceptual viewpoint, Lagrangian averaging provides a powerful tool to
study wave–mean-flow interactions. This is because the material conservation
of key fields (scalar concentrations, vorticity vector, circulation, potential
vorticity) is naturally inherited by the corresponding Lagrangian mean fields.
As a result, the Lagrangian mean of the dynamical equations is often simpler
and more meaningful than the Eulerian mean (Soward, 1972; Bretherton, 1971;
Andrews & McIntyre, 1978; Bühler, 2014; Gilbert & Vanneste, 2018). Relatedly,
the Lagrangian mean emerges naturally in the asymptotic derivation of
wave–averaged models (Grimshaw, 1975; Wagner & Young, 2015). A striking
example of the dynamical relevance of the Lagrangian mean is provided by the
observation that geostrophic balance, the dominant balance in rapidly rotating
flows, continues to hold in the presence of strong waves provided it is
formulated in terms of Lagrangian mean velocity and pressure instead of their
instantaneous or Eulerian mean values (Moore, 1970; Bühler & McIntyre, 1998;
Kafiabad et al., 2021; Kafiabad, 2022).

Despite its advantages, Lagrangian averaging is not widely used as a practical
tool, mainly because Lagrangian means are difficult to compute numerically.
Most numerical models are intrinsically Eulerian and provide the fields of
interest at fixed spatial locations, typically grid points. The standard
approach for the computation of Lagrangian means is then to seed a large
number of passive particles in the flow, track them (forward or backward in
time) using interpolated velocities, and apply time averaging to the resulting
Lagrangian time series (e.g. Nagai et al., 2015; Shakespeare & Hogg, 2017,
2018, 2019; Shakespeare et al., 2021; Bachman et al., 2020; Jones et al.,
2022). This has a high computational cost, requires a large memory allocation,
suffers from possible particle clustering and, as discussed in Kafiabad
(2022), is difficult to parallelise efficiently (see Shakespeare et al. (2021)
for a parallel implementation).

To circumvent the difficulties of particle tracking, Kafiabad (2022) developed
a grid-based method that computes the Lagrangian mean directly on an Eulerian
grid, building the mean through a time step iteration carried out over
successive averaging intervals. By eliminating the need to compute explicit
particle trajectories, the method reduces memory demands and simplifies
integration into parallelised numerical models. The present paper starts with
the recognition that the algorithm of Kafiabad (2022) is a particular
discretisation of a PDE governing the evolution of what we term partial
Lagrangian mean, that is, the mean carried out only up to some intermediate
time in the averaging interval. We formulate this PDE using the position of
particles at the intermediate time as independent spatial variable, as in
Kafiabad (2022). The (total) Lagrangian mean is then obtained by taking the
intermediate time to be the end of the averaging interval.

In this form, the Lagrangian mean does not match Andrews & McIntyre’s (1978)
definition of the generalised Lagrangian mean (GLM): this requires the mean
fields to be expressed as functions of the mean position of fluid particles.
To achieve this, it is necessary to relate the mean positions of particles to
their positions at the end of averaging interval, and to carry out a remapping
of the Lagrangian mean fields. This constitutes our strategy 1 for the
computation of generalised Lagrangian means. We show that the algorithm of
Kafiabad (2022) amounts to a semi-Lagrangian discretisation of the PDEs of
strategy 1. We propose an alternative strategy, strategy 2, which formulates
PDEs directly for the partial Lagrangian means using the mean position as
independent spatial variable. The PDEs involved in both strategies can be
solved by broad classes of numerical methods: finite differences, finite
volumes, finite elements or spectral methods. We illustrate this with a
pseudospectral Fourier implementation for a shallow-water flow in a doubly
periodic domain.

The paper is structured as follows. We introduce notation and define the
Lagrangian means in §2. We derive the PDEs of the two strategies in §3. We
discuss their numerical implementation and present two application to a
shallow-water simulation in §4. The choice of strategy, their advantages and
costs are discussed in §5. Technical aspects including the averaging of
tensorial fields are relegated to appendices.

##  2 Formulation

We consider fluid motion in a two- or three-dimensional Euclidean space. We
denote the flow map by $\bm{\varphi}$, with
$\bm{\varphi}(\bm{a},t)\in\mathbb{R}^{2}$ or $\mathbb{R}^{3}$ the position at
time $t$ of a particle identified by its label $\bm{a}$ (which can be taken as
the position at $t=0$). The flow map and velocity field are related by

$\partial_{t}\bm{\varphi}(\bm{a},t)=\bm{u}(\bm{\varphi}(\bm{a},t),t).$ (1)

Lagrangian averaging is averaging at fixed particle label $\bm{a}$, in
contrast with Eulerian averaging which fixes the spatial position. Both can
involve different types of means: temporal, spatial or – as often used in
theoretical work – ensemble mean. Here we focus on a straigthforward time
average, of the form

$\bar{g}(\tau)=\frac{1}{T}\int_{\tau}^{\tau+T}g(s)\,\mathrm{d}s$ (2)

when applied to a function $g(t)$ that depends on time only. Eq. (2)
introduces the notation $\tau$ for the time at which the averaging is carried
out and $T$ for the averaging period. Usually the middle of averaging interval
$\tau+T/2$ is used as an argument for the mean function, but we prefer to
adopt $\tau$ (the beginning of averaging interval) to simplify the notation in
the upcoming derivations. A simple shift in time switches from one convention
to the other. A weight function could be inserted in the integrand of (2) to
generalise the definition of the average; this would lead to minimal changes
in what follows.

The Lagrangian mean trajectory associated with (2) is represented by the
Lagrangian mean map ${\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}$
defined by

${\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau)\coloneq\frac{1}{T}\int_{\tau}^{\tau+T}\bm{\varphi}(\bm{a},s)\,\mathrm{d}s.$
(3)

Thus ${\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau)$
returns the mean position from $\tau$ to $\tau+T$ of the particle labelled by
$\bm{a}$. The definition (3) makes sense in $\mathbb{R}^{n}$, when
${\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}$
can be interpreted as a vector and averaged component-wise, but not on other
manifolds where more complicated definitions are necessary (Gilbert &
Vanneste, 2018). The (generalised) Lagrangian mean of a scalar function
$f(\bm{x},t)$ is then defined by

${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}({\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau),\tau)\coloneq\frac{1}{T}\int_{\tau}^{\tau+T}f(\bm{\varphi}(\bm{a},s),s)\,\mathrm{d}s.$
(4)

Hence ${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)$
is the average of $f$ along the trajectory of the fluid parcel, regarded as a
function of the Lagrangian mean position $\bm{x}$ and time $\tau$.

Our aim is the development of an efficient numerical method for the
computation of ${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}$
that relies on solving PDEs, which can be discretised in a variety of ways,
rather than on tracking ensembles of particle trajectories. We propose two
strategies and derive the corresponding PDEs in the next section.

##  3 Two strategies

(a)

$t$$f$${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{\varphi}(\bm{a},\tau+T),\tau)$$={\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}({\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau),\tau)$$\tau$$\tau+T$$\tilde{f}(\bm{\varphi}(\bm{a},t),t;\tau)$$f(\bm{\varphi}(\bm{a},t),t)$

(b)

$t$$\bm{x}$$\bm{\varphi}(\bm{a},\tau+T)$${\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau)$$\tau$$\tau+T$$\bm{\varphi}(\bm{a},t)$$\bar{\bm{\varphi}}(\bm{a},t;\tau)$

Figure 2: Lagrangian means and partial Lagrangian means in an averaging
interval $(\tau,\tau+T)$: (a) means of a function $f$ evaluated along the
trajectory $\bm{\varphi}(\bm{a},t)$ of a fluid parcel labelled by $\bm{a}$,
and (b) means of the position, shown here along a single coordinate axis.

Following Kafiabad (2022), we introduce another representation of the
Lagrangian mean via tilde functions defined by

$\displaystyle{\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{\varphi}(\bm{a},\tau+T),\tau)$
$\displaystyle\coloneq\frac{1}{T}\int_{\tau}^{\tau+T}f(\bm{\varphi}(\bm{a},s),s)\,\mathrm{d}s$
(5a) $\displaystyle={\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}({\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau),\tau).$
(5b)

Comparing this definition with (4) shows that the overbar indicates a mean
along a trajectory identified by the Lagrangian mean position, whereas the
tilde indicates the same mean but with the trajectory identified by the actual
position at $t=\tau+T$, that is, at the end of the averaging. This is
illustrated in figure 2. We also introduce the partial mean versions of (3),
(4) and (5b), namely

$\displaystyle\bar{\bm{\varphi}}(\bm{a},t;\tau)$
$\displaystyle\coloneq\frac{1}{t-\tau}\int_{\tau}^{t}\bm{\varphi}(\bm{a},s)\,\mathrm{d}s,$
(6a) $\displaystyle\bar{f}(\bar{\bm{\varphi}}(\bm{a},t;\tau),t;\tau)$
$\displaystyle\coloneq\frac{1}{t-\tau}\int_{\tau}^{t}f(\bm{\varphi}(\bm{a},s),s)\,\mathrm{d}s,$
(6b) $\displaystyle\tilde{f}(\bm{\varphi}(\bm{a},t),t;\tau)$
$\displaystyle\coloneq\frac{1}{t-\tau}\int_{\tau}^{t}f(\bm{\varphi}(\bm{a},s),s)\,\mathrm{d}s,$
(6c)

as in Kafiabad (2022); see figure 2. Clearly the partial means give the total
means when evaluated at $t=\tau+T$. We emphasise that all means used in the
paper are Lagrangian means and that we only indicate this explicitly by a
superscript L for the total means, to distinguish them from the partial means
which are undecorated. The counterpart of (5b) holds for the partial means:

$\tilde{f}(\bm{\varphi}(\bm{a},t),t;\tau)=\bar{f}(\bar{\bm{\varphi}}(\bm{a},t;\tau),t;\tau).$
(7)

Since time average quantities vary over time scales larger than the averaging
period $T$, it is neither necessary nor desirable to compute Lagrangian means
at each of the times at which the velocity $\bm{u}$ and scalar field $f$ are
known, typically discrete times separated by a small time step. Rather, we
think of the averaging time $\tau$ as a slow variable and propose to compute
the Lagrangian means only at $\tau=\tau_{n}=nT$ for $n=0,1,2,\cdots$. We can
therefore carry out independent computations for each $t_{n}$, each involving
only the fields for $t\in(\tau_{n},\tau_{n}+T)$. We now focus on one such
interval and, to lighten the notation, drop the parametric dependence on
$\tau$ from the partial means in (6). For instance, we use
$\bar{f}(\bar{\bm{\varphi}}(\bm{a},t),t)$ instead of
$\bar{f}(\bar{\bm{\varphi}}(\bm{a},t;\tau),t;\tau)$ in the following
derivations, keeping in mind the now implicit dependence of $\bar{f}$ and
$\bar{\bm{\varphi}}$ on $\tau$. A warning about notation might be useful: in
what follows, we use the symbol $\bm{x}$ as a generic dummy variable, without
attributing it the specific meaning of either an actual or Lagrangian mean
position. This meaning is determined by the function in which $\bm{x}$
appears. Thus $\bm{x}$ in $\tilde{f}(\bm{x},t)$ is interpreted as an actual
position at time $t$, whereas in $\bar{f}(\bm{x},t)$ it is interpreted as a
(partial) Lagrangian mean position.

We now formulate two distinct strategies for the computation of the Lagrangian
mean ${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)$.

###  3.1 Strategy 1

Our first strategy parallels that of Kafiabad (2022) and consists in solving a
PDE for $\tilde{f}(\bm{x},t)$, evaluating the result at $t=\tau+T$ to obtain
${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)$, then
deducing ${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)$
by applying a suitable re-mapping. To derive the PDE for $\tilde{f}(\bm{x},t)$
we take the time derivative of (6c) at fixed label $\bm{a}$ and use the chain
rule and (1) to find

$\displaystyle\partial_{t}\tilde{f}(\bm{\varphi}(\bm{a},t),t)+\bm{u}(\bm{\varphi}(\bm{a},t),t)\bcdot\bnabla\tilde{f}(\bm{\varphi}(\bm{a},t),t)=\frac{1}{t-\tau}\left(f(\bm{\varphi}(\bm{a},t),t)-\tilde{f}(\bm{\varphi}(\bm{a},t),t)\right),$
(8)

Here and in what follows, the gradient $\bnabla$ is taken with respect to the
first argument of $\tilde{f}$. Replacing $\bm{\varphi}(\bm{a},t)$ by $\bm{x}$
as independent variable yields the sought PDE,

$\partial_{t}\tilde{f}(\bm{x},t)+\bm{u}(\bm{x},t)\bcdot\bnabla\tilde{f}(\bm{x},t)=\frac{f(\bm{x},t)-\tilde{f}(\bm{x},t)}{t-\tau},$
(9)

which can be integrated in from $\tau$ to $t=\tau+T$ to find the total mean
${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)=\tilde{f}(\bm{x},\tau+T;\tau)$.
This is a forced advection equation in which the forcing can be interpreted as
a time-dependent relaxation of $\tilde{f}$ to $f$. In a bounded domain, the
solution of (9) requires no boundary conditions since the differentiation
$\bm{u}(\bm{x},t)\bcdot\bnabla$ is along the boundary. The initial condition
is that $\tilde{f}(\bm{x},\tau)=f(\bm{x},\tau)$ so that the right-hand side is
finite.

label spacephysical space
$\subseteq\mathbb{R}^{n}$$\bm{\varphi}$$\bm{a}$$\bar{\bm{\varphi}}(\bm{a},t;\tau)$$\bm{x}$$\bar{\bm{\varphi}}$${\bm{\varXi}}$${\bm{\varXi}^{-1}}$
Figure 3: Flow map $\bm{\varphi}$, Lagrangian (partial) mean map
$\bar{\bm{\varphi}}$, lift map $\bm{\varXi}$ and its inverse
$\bm{\varXi}^{-1}$.

Computing ${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}$ may be all that is
needed for applications in which the spatial distribution of the Lagrangian
mean is not important. For example, wave-averaged geostrophy – the modified
form of geostrophic balance expressed in terms of Lagrangian mean quantities –
can be validated by comparing Lagrangian mean velocity and pressure gradient
at the same position (equivalently, the same label) regardless of where the
position is located in physical space (Kafiabad et al., 2021; Kafiabad, 2022).

However, in many other applications, it is necessary to compute the
(generalised) Lagrangian mean field ${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},t)$
for specified Lagrangian mean positions $\bm{x}$. This (generalised)
Lagrangian mean field can be deduced from ${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},t)$ by
considering the lift map $\bm{\varXi}(\bm{x},t)$ which returns the position at
time $t$ of the particle with $\bm{x}$ as partial mean position from $\bar{t}$
to $t$, that is,

$\quad\bm{\varXi}(\bar{\bm{\varphi}}(\bm{a},t),t)=\bm{\varphi}(\bm{a},t)$ (10)

(Andrews & McIntyre, 1978; Bühler, 2014). Its inverse
$\bm{\varXi}^{-1}(\bm{x},t)$ returns the partial mean position of the particle
that passes through $\bm{x}$ at $t$:

$\bm{\varXi}^{-1}(\bm{\varphi}(\bm{a},t),t)=\bar{\bm{\varphi}}(\bm{a},t)=\frac{1}{t-\tau}\int_{\tau}^{t}\bm{\varphi}(\bm{a},s)\,\mathrm{d}s,$
(11)

where we used the definition (6a) of the mean position. The map $\bm{\varXi}$
and its inverse $\bm{\varXi}^{-1}$ are depicted in the figure 3. The relation
(5b) between ${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}$ and
${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}$
can then be written in terms of $\bm{\varXi}^{-1}$ as

${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)={\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{\varXi}^{-1}(\bm{x},\tau+T),\tau).$
(12)

Now, comparing (11) with (6c) makes it clear that the components of
$\bm{\varXi}^{-1}$ can be viewed as instances of functions $\tilde{f}$ with
$f(\bm{x})=x_{i}$. Thus, we can rewrite (9) for this special case to obtain

$\partial_{t}\bm{\varXi}^{-1}(\bm{x},t)+\bm{u}(\bm{x},t)\bcdot\bnabla\bm{\varXi}^{-1}(\bm{x},t)=\frac{\bm{x}-\bm{\varXi}^{-1}(\bm{x},t)}{t-\tau}.$
(13)

Alternatively, we can take the time derivative of (11) and replace
$\bm{\varphi}(\bm{a},t)$ by $\bm{x}$ to arrive at (13). Integrating (13)
provides the means to effect the remapping (12) between ${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}$ and
${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}$.

To recapitulate, our first strategy consists in solving the PDEs (9) and (13)
from $\tau$ to $\tau+T$ to obtain ${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)=\tilde{f}(\bm{x},\tau+T;\tau)$
and $\bm{\varXi}(\bm{x},\tau+T;\tau)$, then using (12) to compute
${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}$
by interpolation. The algorithm proposed by Kafiabad (2022) turns out to be a
particular discretisation of this strategy (see §4.1 below).

###  3.2 Strategy 2

Our second strategy bypasses the use of $\tilde{f}$ and is instead based on
PDEs for $\bar{f}$ and $\bm{\varXi}$. To derive these we first note that
taking the time derivative of (6a) gives

$\partial_{t}\bar{\bm{\varphi}}(\bm{a},t)=\frac{\bm{\varphi}(\bm{a},t)-{\bar{\bm{\varphi}}(\bm{a},t)}}{t-\tau}\eqcolon\bar{\bm{u}}(\bar{\bm{\varphi}}(\bm{a},t),t).$
(14)

The second equality defines the auxiliary velocity field $\bar{\bm{u}}$ as the
time derivative of the partial Lagrangian mean position. Using (10), this
velocity field can be written in terms of the lift map as

$\bar{\bm{u}}(\bm{x},t)=\frac{\bm{\varXi}(\bm{x},t)-\bm{x}}{t-\tau},$ (15)

where the dummy variable $\bm{x}$ can be thought of as the partial mean
position. We emphasise that $\bar{\bm{u}}(\cdot,t)$, like
$\bm{\varXi}(\cdot,t)$, depends implicitly on $\tau$ and warn that it should
not be interpreted as the partial mean of the Lagrangian velocity: as
discussed in appendix A, its value at the end of the averaging interval, for
$t=\tau+T$, differs from the usual Lagrangian mean velocity, that is, the time
derivative of ${\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau)$
with respect to $\tau$.

Now, differentiating (6b) with respect to $t$ at fixed label $\bm{a}$ and
using (14) leads to

$\partial_{t}\bar{f}(\bar{\bm{\varphi}}(\bm{a},t),t)+\bar{\bm{u}}(\bar{\bm{\varphi}}(\bm{a},t),t)\bcdot\bnabla\bar{f}(\bar{\bm{\varphi}}(\bm{a},t),t)=\frac{f(\bm{\varphi}(\bm{a},t),t)-\bar{f}(\bar{\bm{\varphi}}(\bm{a},t),t)}{t-\tau}.$
(16)

We obtain the desired PDE for $\bar{f}(\bm{x},t)$ by using (10) and replacing
$\bar{\bm{\varphi}}(\bm{a},t)$ by the independent variable $\bm{x}$ to write

$\partial_{t}\bar{f}(\bm{x},t)+\bar{\bm{u}}(\bm{x},t)\bcdot\bnabla\bar{f}(\bm{x},t)=\frac{f(\bm{\varXi}(\bm{x},t),t)-\bar{f}(\bm{x},t)}{t-\tau}.$
(17)

This is a forced advection equation, analogous to the PDE (9) governing
$\tilde{f}(\bm{x},t)$. However, unlike (9) it is not closed since it involves
$\bm{\varXi}(\bm{x},t)$, explicitly on the right-hand side and implicitly
through $\bar{\bm{u}}(\bm{x},t)$ on the left-hand side. It needs to be solved
along an equation for $\bm{\varXi}(\bm{x},t)$. We derive this equation by
taking the time derivative of (10) at fixed $\bm{a}$ and using (14) and (1) to
obtain

$\partial_{t}\bm{\varXi}(\bar{\bm{\varphi}}(\bm{a},t),t)+\bar{\bm{u}}(\bar{\bm{\varphi}}(\bm{a},t),t)\bcdot\bnabla\bm{\varXi}(\bar{\bm{\varphi}}(\bm{a},t),t)=\bm{u}(\bm{\varphi}(\bm{a},t),t).$
(18)

Hence, replacing $\bar{\bm{\varphi}}(\bm{a},t)$ by $\bm{x}$ and using (10),

$\partial_{t}\bm{\varXi}(\bm{x},t)+\bar{\bm{u}}(\bm{x},t)\bcdot\bnabla\bm{\varXi}(\bm{x},t)=\bm{u}(\bm{\varXi}(\bm{x},t),t).$
(19)

Strategy 2 consists in solving (17) and (19), with $\bar{\bm{u}}(\bm{x},t)$
defined by (15), for $t\in(\tau,\tau+T)$, then deduce the Lagrangian mean of
$f$ as ${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)=\bar{f}(\bm{x},\tau+T;\tau)$.

The initial conditions for (17) and (19) are that
$\bar{f}(\bm{x},\tau)=f(\bm{x},\tau)$ and $\bm{\varXi}(\bm{x},\tau)=\bm{x}$.
The boundary conditions are non-trivial: in bounded domains,
$\bar{f}(\bm{x},t)$ and $\bm{\varXi}(\bm{x},t)$ are defined on the image of
the label space by the Lagrangian mean map $\bar{\bm{\varphi}}$ (equivalently,
the image of the fluid domain by $\bm{\varXi}^{-1}$). Thus, the problem in
principle involves a boundary moving with velocity $\bar{\bm{u}}$ and can
therefore be difficult to discretise. The common situation where the physical
domain has boundaries that coincide with constant-coordinate surfaces (curves
in two dimensions) is straightforward, however, because the component-wise
definition of $\bar{\bm{\varphi}}$ in (3) ensures that it maps such boundaries
to themselves so the domain remains fixed. The case of periodic domains is
also straightforward.

##  4 Numerical implementation

The set of equations for each strategy of the previous section can be
discretised in a variety of ways. Here we focus on a pseudospectral
discretisation which we apply to the computation of Lagrangian means in a
turbulent shallow-water flow interacting with a Poincaré wave. We make general
remarks about the choice of strategy and numerical discretisation but leave a
more complete analysis of numerical error and convergence for future studies.

###  4.1 Strategy 1

To solve (9) it is convenient to introduce
$\tilde{g}(\bm{x},t)=(t-\tau)\tilde{f}(\bm{x},t)$, leading to

$\left(\partial_{t}+\bm{u}(\bm{x},t)\bcdot\bnabla\right)\tilde{g}(\bm{x},t)=f(\bm{x},t).$
(20)

Integrating this equation over the averaging period then yields the Lagrangian
mean ${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)=\tilde{g}(\bm{x},\tau+T)/T$.
For simple geometries, periodic in particular, standard pseudospectral methods
provide efficient solvers for (20) and, if the remapping to Lagrangian mean
positions is desired, (13). This is particularly convenient if the fluid’s
governing equations are also solved pseudospectrally, because $f$ is then
available on spectral grid points to evaluate the right-hand sides of (20) and
(13), and on physical grid points to evaluate the nonlinear terms. An
alternative is a semi-Lagrangian discretisation, which leads to the algorithm
of Kafiabad (2022) as detailed in appendix B.

Figure 4: Vorticity field and its Lagrangian mean for a two-dimensional
incompressible inviscid flow: (a) instantaneous vorticity $\zeta$ at $t=25$,
and (b–d) $\tilde{\zeta}^{\scriptscriptstyle\mathrm{L}}$ computed by time-
integration of (9) from $\tau=0$ to $T=25$ using (b) a semi-Lagrangian
discretisation with linear interpolation, (c) semi-Lagrangian discretisation
with cubic interpolation and (d) pseudospectral discretisation.

To investigate the validity of numerical solutions of (9) we consider a two-
dimensional incompressible inviscid flow for which the vorticity, $\zeta$ say,
is conserved materially. This implies that

$\zeta(\bm{x},\tau+T)={\tilde{\zeta}\mkern
1.5mu}{\vphantom{\bar{\zeta}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau).$
(21)

Hence, we can calculate ${\tilde{\zeta}\mkern
1.5mu}{\vphantom{\bar{\zeta}}}^{\scriptscriptstyle\mathrm{L}}$ by integrating
(9) (or (20)) from $\tau$ to $\tau+T$ and compare it with the instantaneous
vorticity at $\tau+T$ to study the accuracy of the computed Lagrangian mean.
Note that this is simply a test for (9) using the material conservation of
$\zeta$ as opposed to an application of the Lagrangian mean. As mentioned
earlier, (9) should usually by solved in tandem with (13) to get a meaningful
spatial distribution of Lagrangian mean quantities.

We perform a numerical simulation of this two-dimensional flow without
viscosity in a doubly periodic, $\left[0,\ 2\pi\right]^{2}$ domain, using a
standard pseudospectral discretisation and $2/3$ de-aliasing with $128^{2}$
grid points. We start the simulation with the vorticity

$\zeta(x,y,t=0)=\mathrm{e}^{-(x-\pi+0.1)^{2}-(y-\pi+\pi/3)^{2}}+\mathrm{e}^{-(x-\pi-0.1)^{2}-(y-\pi-\pi/3)^{2}},$
(22)

corresponding to two like-signed vortices which subsequently merge. We use
Heun’s method for the time integration of the governing vorticity equation and
for (20), with time step $\Delta t=0.005$. Figure 4 displays the instantaneous
vorticity $\zeta$ at $t=25$ and ${\tilde{\zeta}\mkern
1.5mu}{\vphantom{\bar{\zeta}}}^{\scriptscriptstyle\mathrm{L}}$ obtained for
$\tau=0$ and $T=25$. As expected from (21), the Lagrangian mean
${\tilde{\zeta}\mkern
1.5mu}{\vphantom{\bar{\zeta}}}^{\scriptscriptstyle\mathrm{L}}$ matches the
instantaneous vorticity $\zeta$. The pseudospectral solution for
${\tilde{\zeta}\mkern
1.5mu}{\vphantom{\bar{\zeta}}}^{\scriptscriptstyle\mathrm{L}}$, shown in panel
(d), in particular, shows an excellent agreement with $\zeta$ in panel (a).
The results of the semi-Lagrangian algorithm of Kafiabad (2022) with,
respectively, linear and cubic interpolations, are shown in panels (b) and (c)
(see appendix B). These show a poorer agreement with panel (a), especially
with the linear interpolation, because of an accumulation of interpolation
errors. The computation reported in figure 4 is however rather extreme in both
the coarseness of the resolution and the length of the averaging interval. We
have confirmed that the three numerical solutions for ${\tilde{\zeta}\mkern
1.5mu}{\vphantom{\bar{\zeta}}}^{\scriptscriptstyle\mathrm{L}}$ converge to
each other and to $\zeta$ as the spatial resolution increases or the length of
the averaging interval decreases (the interested reader will find a
demonstration of this convergence in the supplementary material). Since we
solve the dynamical and Lagrangian-mean equations without explicit
dissipation, both require small time steps. Our investigation for this
particular example shows that the Lagrangian mean equations are less
restrictive than the dynamical equations for the size of time step.

The pseudospectral method leads to the more the accurate results, but it is
not as stable as its semi-Lagrangian counterpart and therefore requires
smaller time steps. The difference arises because the implicit time
integration and numerical smoothing due to interpolation that are inherent to
semi-Lagrangian methods have a stabilising effect.

###  4.2 Strategy 2

We now implement strategy 2 which uses the Lagrangian mean position $\bm{x}$
as independent spatial variable. In the periodic domain we consider, there are
no difficulties associated with moving boundaries and a pseudospectral
discretisation is straightforward. It is convenient to rewrite the PDEs to be
integrated, (17) and (19), in terms of the displacement map

$\bm{\xi}(\bm{x},t)=\bm{\varXi}(\bm{x},t)-\bm{x},$ (23)

since $\bm{\xi}$ is periodic, unlike $\bm{\varXi}$. (This is the partial-mean
analogue of the displacement introduced by Andrews & McIntyre (1978).) As in
strategy 1, it is also convenient to solve for
$\bar{g}(\bm{x},t)=(t-\tau)\bar{f}(\bm{x},t)$ instead of $\bar{f}$. With these
transformations, $\bar{\bm{u}}=\bm{\xi}/(t-\tau)$ and (17) and (19) are
rewritten as

$\displaystyle\partial_{t}\bm{\xi}(\bm{x},t)+\frac{\bm{\xi}(\bm{x},t)}{t-\tau}\bcdot\bnabla\bm{\xi}(\bm{x},t)$
$\displaystyle=\bm{u}(\bm{x}+\bm{\xi}(\bm{x},t),t)-\frac{\bm{\xi}(\bm{x},t)}{t-\tau},$
(24a)
$\displaystyle\partial_{t}\bar{g}(\bm{x},t)+\frac{\bm{\xi}(\bm{x},t)}{t-\tau}\bcdot\bnabla\bar{g}(\bm{x},t)$
$\displaystyle=f(\bm{x}+\bm{\xi}(\bm{x},t),t).$ (24b)

These are the PDEs we solve numerically. When discretising in time, we found
it beneficial for stability to first update $\bm{\xi}$ using (24a), then use
the updated $\bm{\xi}$ for the time integration of (24b).

Below we apply strategy 2 to two examples of rotating shallow-water flows. We
write the governing equations in a non-dimensional form, using a
characteristic length $L$, characteristic velocity $U$, time $L/U$ and mean
height $H$ for scaling, leading to

$\displaystyle\frac{\partial\bm{u}}{\partial
t}+\bm{u}\bcdot\bnabla\bm{u}+\frac{1}{\rm{Ro}}\hat{\bm{z}}\times\bm{u}$
$\displaystyle=-\frac{1}{\rm{Fr}^{2}}\ \bnabla
h+\frac{1}{\rm{Re}}\nabla^{2}\bm{u},$ (25a) $\displaystyle\frac{\partial
h}{\partial t}+\bnabla\bcdot(h\bm{u})$ $\displaystyle=0\ ,$ (25b)

where we introduce the standard dimensionless numbers (e.g. Vallis, 2017)

${\rm{Ro}}=\frac{U}{fL}\
,\quad{\rm{Fr}}=\frac{U}{\sqrt{gH}}\quad\textrm{and}\quad{\rm{Re}}=\frac{UL}{\nu}.$
(26)

In the above, $g$ is the gravitational acceleration, $\nu$ the kinematic
viscosity, $f$ the Coriolis parameter and $\hat{\bm{z}}$ the vertical unit
vector.

####  4.2.1 Interaction of a turbulent geostrophic flow with a strong wave

Figure 5: Initial condition of the shallow-water simulation: (a) vertical
vorticity of the geostrophic flow and (b) height field of the Poincaré wave.

In the first example, we compute Lagrangian means in a simulation of a
turbulent flow interacting with a Poincaré wave in a rotating shallow water.
We initialise our simulation with a turbulent flow that is initially in
geostrophic balance, with vorticity $\zeta=\partial_{x}v-\partial_{x}u$ shown
in figure 5(a), and superimpose a mode-1 Poincaré wave, with the height field
shown in figure 5(b). The initial geostrophic flow is generated from the
output of an incompressible two-dimensional Navier-Stokes simulation, which
has reached a fully-developed turbulent state, and the height field that is in
geostrophic balance with this vortical flow. We use the root-mean square
velocity of the geostrophic flow as the characteristic velocity $U$, and
choose the length scale of first Fourier mode as characteristic length $L$.
This makes the dimensionless doubly periodic domain $\left[0,\
2\pi\right]^{2}$ which we discretise with $256\times 256$ grid points. The
right-travelling mode-1 Poincaré wave has the form

$u^{\prime}=a\cos(x-\omega t),\quad
v^{\prime}=\frac{a}{\omega\rm{Ro}}\sin(x-\omega t),\quad
h^{\prime}=\frac{a}{\omega}\cos(x-\omega t),$ (27)

with the intrinsic frequency
$\omega=\left(\rm{Ro}^{-2}+\rm{Fr}^{-2}\right)^{1/2}$ and the velocity
amplitude $a$ is a constant, taken as $a=-1.8$ in our simulation. This implies
wave velocities that are almost twice as large as the geostrophic velocities.
It is in this sense that we regard the wave as strong. We set the
dimensionless parameters to

${\rm{Ro}}=0.1\ ,\quad{\rm{Fr}}=0.5\quad\textrm{and}\quad{\rm{Re}}=3.84\times
10^{3},$ (28)

which results in $\omega=10.2$.

We evaluate (27) at $t=0$ and add the wave fields to the geostrophic field to
form the initial condition. We solve the dynamical equations (25) in tandem
with the Lagrangian mean equations (24) over a single averaging time interval
taken to be $T=2.2$, corresponding to approximately $3.6$ wave periods. We use
a de-aliased pseudospectral discretisation and a forward Euler integrator,
with the time step of $1.25\times 10^{-4}$ for (25) and $2.5\times 10^{-4}$
for (24). A bilinear interpolation is used to evaluate $\bm{u}$ and $f$ at
$\bm{x}+\bm{\xi}(\bm{x},t)$ in the right-hand sides of (24). The link for the
scripts and data used to produce the results of this section is provided in
the “Data availability statement” at the end of this paper.

Figure 6: Potential vorticity $q$ and its Lagrangian and Eulerian means in the
shallow-water simulation: (a) instantaneous PV $q$ at $t=1.1$, (b) Lagrangian
mean $\bar{q}^{\scriptscriptstyle\mathrm{L}}$, and (c) Eulerian mean PV. The
mean fields are averaged from $\tau=0$ to $T=2.2$. All panels share the same
colourbar.

Figure 1, briefly discussed in §1, displays the vorticity field $\zeta$ at
$t=T/2$ (panel (a)) and its Lagrangian and Eulerian means (panels (b) and
(c)). The strong wave dominates the instantaneous vorticity field. It is
filtered successfully by time averaging. Clearly, the Lagrangian mean captures
small-scale structures in the vorticity which are blurred by the Eulerian
mean. This blurring is the result of the advection of the vorticity by the
velocity field associated with the wave which causes the vorticity structures
to oscillate with the wave period. By construction, the Lagrangian mean
removes these oscillations, leading to a sharper definition of the flow
features.

It is interesting to examine the effect of Lagrangian and Eulerian averaging
on the potential vorticity (PV) $q=(1+\rm{Ro}\ \zeta)/h$. In the absence of
dissipation ($\rm{Re}\to\infty$), this is a materially conserved field,
meaning that $q(\bm{\varphi}(\bm{a},t),t)=q_{0}(\bm{a})$, with $q_{0}$
determined by the initial condition. By definition (4), the Lagrangian mean PV
then satisfies
$\bar{q}^{\scriptscriptstyle\mathrm{L}}(\bar{\bm{\varphi}}(\bm{a},t),t)=q_{0}(\bm{a})$.
Thus both $q$ and $\bar{q}^{\scriptscriptstyle\mathrm{L}}$ are (smooth) re-
arrangements of the initial PV and hence re-arrangements of one another,
specifically
$\bar{q}^{\scriptscriptstyle\mathrm{L}}(\bm{x},t)=q(\bm{\varXi}(\bm{x},t),t)$.
This imposes constraints such as the two fields sharing the same values for
their local extrema. Because $\bar{\bm{\varphi}}$ is not area-preserving, the
distribution functions of $q$ and $\bar{q}^{\scriptscriptstyle\mathrm{L}}$
(measuring the area of regions where the fields are below specified values) do
not coincide. Figure 6 shows $q$ at $t=T/2$ and
$\bar{q}^{\scriptscriptstyle\mathrm{L}}$ as well as the Eulerian mean PV. The
Lagrangian mean PV $\bar{q}^{\scriptscriptstyle\mathrm{L}}$ appears as a
slight deformation of $q$, consistent with it being rearrangement by a map
$\bm{\varXi}$ that is close to the identity. In contrast, the Eulerian mean
PV, which is not materially transported, shows blurred features, with in
particular extrema that are substantially reduced compared with those of $q$
and $\bar{q}^{\scriptscriptstyle\mathrm{L}}$. There is a strong argument that
the study of wave–mean-flow interactions, in the shallow-water model and more
broadly, would benefit from the systematic analysis of Lagrangian mean fields
such as the ones displayed in panels (b) of figures 1 and 6.

####  4.2.2 Moderate-Rossby-number flow initialised in geostrophic balance

Figure 7: (a) Instantaneous vorticity $\zeta(\bm{x},t)$ at $t=1.1$; (b)
Lagrangian mean vorticity
$\bar{\zeta}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau=0)$ corresponding to
averaging from $t=0$ to $T=2.2$. (c) Eulerian mean vorticity averaged from
$t=0$ to $2.2$; (d) Enstrophy spectra of instantaneous velocity at $t=2.2$ and
the Lagrangian and Eulerian mean fields; (e) Lagrangian perturbation field;
(f) Eulerian perturbation field.

In the second example, we initialise the shallow-water simulation with the
same geostrophic flow as in §4.2.1 but without any waves. The set-up is
identical, but with a resolution of $128^{2}$ and an increased Rossby number
$\rm{Ro}=0.6$. As the flow evolves at this relatively high $\rm{Ro}$, small
scale imbalance is generated as shown in the instantaneous vorticity of panel
(a) in figure 7 and the associated enstrophy spectrum in panel (d). Panels (b)
and (e) show the Lagrangian mean vorticity field and the corresponding
perturbation field, respectively. The perturbation field is computed as
${\mathrel{\mathop{\zeta}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\zeta}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)-\zeta(\bm{x}+\bm{\xi}(\bm{x},\tau+T),\tau+T)$
with, in this case, $\tau=0$ and $T=2.2$. The Lagrangian average preserves the
large-scale balanced flow so that the corresponding perturbation extracts the
small-scale imbalance. This is evident in the enstrophy spectra of panel (d),
where the instantaneous and Lagrangian mean spectra overlap at small
wavenumbers and the jump in the tail of the instantaneous spectrum is removed
in the Lagrangian mean spectrum. In contrast, the mean flow is weakened in the
Eulerian average (panel (c)) because, at this moderate $\rm{Ro}$, the vortices
are advected substantially during the averaging period. This is corroborated
in the Eulerian mean enstrophy spectrum, which has lower values for small
wavenumbers, and follows the jump of the instantaneous spectrum at large
wavenumbers. Consequently, balanced flow features dominate the Eulerian
perturbation field (panel (f)). The results of this figure are similar to
those of the synthetic flow considered by Shakespeare et al. (2021).

##  5 Discussion

This paper presents a novel approach for the numerical computation of
Lagrangian means which relies on solving PDEs rather than tracking particles.
We propose two strategies, each leading to a separate set of PDEs. Both
strategies are based on the derivation of equations governing the evolution of
partial means with respect to $t$. These partial means are defined as averages
over a subset $(\tau,t)$ of each averaging interval $(\tau,\tau+T)$ and yield
the (total) means for $t=\tau+T$. Strategy 1 uses the position of particles at
time $t$ as independent spatial variable. Hence, it requires a map from the
positions at the final time $t=\tau+T$ to the Lagrangian mean positions to
ultimately present the results in terms of the latter, as is standard in GLM
theory. Strategy 2 directly computes the Lagrangian means in terms of
Lagrangian mean positions.

A natural question is which of the two strategies should be preferred. There
is no definitive answer: each strategy has pros and cons. If the spatial
distribution of Lagrangian means does not matter in an application, it
suffices to solve Eq. (9) of strategy 1. When the spatial distribution is
needed, strategy 1 requires the re-mapping (12) which can be affected by
clustering: the mean positions $\bm{\varXi}^{-1}(\bm{x},\tau+T)$ obtained from
(13) for $\bm{x}$ on a regular grid may have a highly non-uniform
distribution. This can lead to large numerical errors in the interpolation
required to discretise (12). Strategy 2 circumvents this issue, as the
generalised Lagrangian mean is computed directly on the desired spatial grid
points. However, this advantage comes at the computational cost of evaluating
more complicated right-hand sides in (17) and (19), which require
interpolation at each averaging time step. Furthemore, strategy 2 leads to
PDEs posed on a moving domain, unless the domain is periodic or has boundaries
that correspond to constant coordinates.

As discussed in Kafiabad (2022), the full potential of our approach in saving
memory and reducing computational cost is realised when the PDEs for the
Lagrangian mean fields are solved on-the-fly, together with the dynamical
model (as opposed to offline, using saved model outputs). In this case, it is
beneficial to solve the Lagrangian mean PDEs using a numerical scheme that
closely matches that of the dynamical model, because the instantaneous values
of $f$ and $\bm{u}$ (required to solve the Lagrangian mean PDEs) are readily
available at the same (physical or spectral) grid points. Moreover, the
reasons that led to a particular choice of numerical discretisation for the
dynamical equations – such as the type of boundary conditions – typically also
apply to the Lagrangian mean PDEs.

In the main body of the paper, we restrict our attention to the Lagrangian
averaging of a scalar function $f(\bm{x},t)$. The averaging of vectors,
differential forms and more general tensors is however of interest in
applications. For instance, the Lagrangian mean of the momentum 1-form
$\bm{u}\bcdot\mathrm{d}\bm{x}$ (the integrand in Kelvin’s circulation) and of
the magnetic flux 2-form play crucial roles in the theory of wave–mean-flow
interactions in fluid dynamics and MHD (Soward, 1972; Andrews & McIntyre,
1978; Holm, 2002; Gilbert & Vanneste, 2018, 2021). The derivations in §3
generalise straightforwardly to tensors when the language of push-forwards,
pull-backs and Lie derivatives is employed. We illustrate this in appendix C
by generalising Eq. (17) of strategy 2 for the partial Lagrangian mean of
$f(\bm{x},t)$ to a tensor field $\tau(\bm{x},t)$.

We conclude by noting that practical averages such as the time average in
(2)–(4) do not satisfy exactly the axioms of the more abstract averages
assumed in the development of GLM and similar theories. In particular, the
basic requirement that averaging leaves mean quantities unchanged, that is,
$\bar{\bar{g}}\not=\bar{g}$, fails for (2), though the difference is small if
there is a clear time-scale separation between mean flow and perturbation.
Interpreting numerically computed Lagrangian mean fields in light of these
theories will therefore require to understand how theoretical predictions are
affected by the precise nature of the average.

##  Appendix A Partial and total Lagrangian mean velocities

We show that the partial mean velocity $\bar{\bm{u}}$ in (15) is not related
to the Lagrangian-mean velocity ${\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}$
as might be expected naively: ${\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)\not=\bar{\bm{u}}(\bm{x},\tau+T;\tau)$,
where we have reinstated the parameteric dependence on $\tau$ for clarity. The
definition of ${\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}$
itself is not without ambiguity: the most natural definition is as the
derivative of the total Lagrangian mean map with respect to the slow time,
i.e.

$\partial_{\tau}{\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau)={\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}({\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau),\tau).$
(29)

Taking the derivative of the definition (3) of
${\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}$
gives the explicit form

${\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}({\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{a},\tau),\tau)=\frac{1}{T}\left(\bm{\varphi}(\bm{a},\tau+T)-\bm{\varphi}(\bm{a},\tau)\right).$
(30)

In contrast, the partial mean velocity $\bar{\bm{u}}$ is defined via a
derivative with respect to the fast time $t$ as seen from the definition (14)
which takes the form

$\partial_{t}\bar{\bm{\varphi}}(\bm{a},t;\tau)\eqcolon\bar{\bm{u}}(\bar{\bm{\varphi}}(\bm{a},t),t;\tau).$
(31)

on reinstating the parametric dependence on $\tau$.

The difference between $\bar{\bm{u}}(\bm{x},\tau+T;\tau)$ and
${\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)$
should not come as a surprise since $\bar{\bm{u}}(\bm{x},\tau+T;\tau)$ is
constructed independently for each value of $\tau$ and hence cannot capture
the change of ${\mathrel{\mathop{\bm{\varphi}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{\varphi}}}}^{\scriptscriptstyle\mathrm{L}}(\cdot,\tau)$
measured by ${\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}$.
It can be made explicit by deducing from (30) that

$\displaystyle{\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)$
$\displaystyle=\frac{1}{T}\left(\bm{\varphi}(\bar{\bm{\varphi}}^{-1}(\bm{x},\tau+T;\tau),\tau+T)-\bm{\varphi}(\bar{\bm{\varphi}}^{-1}(\bm{x},\tau+T;\tau),\tau)\right)$
$\displaystyle=\frac{1}{T}\left(\bm{\varXi}(\bm{x},t;\tau)-\bm{\varphi}(\bar{\bm{\varphi}}^{-1}(\bm{x},\tau+T;\tau),\tau)\right),$
(32)

which differs from
$\bar{\bm{u}}(\bm{x},\tau+T;\tau)=\left(\bm{\varXi}(\bm{x},\tau+T;\tau)-\bm{x}\right)/T$
since
$\bm{\varphi}(\bar{\bm{\varphi}}^{-1}(\bm{x},\tau+T;\tau),\tau)\not=\bm{x}$.

Note that in GLM, ${\mathrel{\mathop{\bm{u}}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\bm{u}}}}^{\scriptscriptstyle\mathrm{L}}$
is defined as the Lagrangian mean of the components of $\bm{u}$ instead of via
(29). In Euclidean space the two definitions are equivalent since

$\frac{1}{T}\int_{\tau}^{\tau+T}\bm{u}(\bm{\varphi}(\bm{a},s),s)\,\mathrm{d}s=\frac{1}{T}\int_{\tau}^{\tau+T}\partial_{s}\bm{\varphi}(\bm{a},s)\,\mathrm{d}s=\frac{\bm{\varphi}(\bm{a},\tau+T)-\bm{\varphi}(\bm{a},\tau)}{T},$
(33)

using (1). This equivalence does not hold for definitions of the Lagrangian
mean flow that have been proposed as geometric alternatives to GLM valid on
any manifold (Soward & Roberts, 2010; Gilbert & Vanneste, 2018; Vanneste &
Young, 2022).

##  Appendix B Semi-Lagrangian discretisation of strategy 1

We show that the algorithm developed by Kafiabad (2022) for the ‘grid-based
calculation of Lagrangian average’ is equivalent to a semi-Lagrangian
discretisation of (9) and (13). We denote the grid points by $\bm{X}_{i}$,
with $i$ a multi-index in dimension more than one, and the set of all grid
points by $\\{\bm{X}_{\alpha}\\}$. The semi-Lagrangian discretisation (9)
gives the value $\tilde{f}(\bm{X}_{i},t_{n+1})$ of $\tilde{f}$ at each grid
point at $t_{n+1}$, the end of a time step, in terms of the set of values
$\tilde{f}(\\{\bm{X}_{\alpha}\\},t_{n})$ at $t_{n}$, the start of the time
step. It is based on a semi-Lagrangian discretisation of the equivalent PDE
(20), namely

$\tilde{g}(\bm{X}_{i},t_{n+1})=\tilde{g}(\bm{Y}_{i},t_{n})+(t_{n+1}-t_{n})f(\bm{X}_{i},t_{n+1}),$
(34)

where we use a backward Euler step for the time stepping. Here $\bm{Y}_{i}$
denotes the position at $t_{n}$ of the particle that passes through
$\bm{X}_{i}$ at $t_{n+1}$; it can be evaluated using

$\bm{Y}_{i}=\bm{X}_{i}-(t_{n+1}-t_{n})\bm{u}(\bm{X}_{i},t_{n+1}),$ (35)

which follows from approximating the derivative in (1) by a backward finite
difference. Rewriting (34) in terms of $\tilde{f}$ and rearranging we obtain

$\tilde{f}(\bm{X}_{i},t_{n+1})=\frac{(t_{n}-\tau)\tilde{f}(\bm{Y}_{i},t_{n})+(t_{n+1}-t_{n})f(\bm{X}_{i},t_{n+1})}{t_{n+1}-\tau},$
(36)

which is the same as Eq. (2.6) in Kafiabad (2022) when the last term is
replaced by their (2.7). Using the trapezoidal rule instead of backward Euler
leads to

$\tilde{f}(\bm{X}_{i},t_{n+1})=\frac{(t_{n}-\tau)\tilde{f}(\bm{Y}_{i},t_{n})+(t_{n+1}-t_{n})\left[f(\bm{Y}_{i},t_{n})+f(\bm{X}_{i},t_{n+1})\right]/2}{t_{n+1}-\tau},$
(37)

which is Eq. (2.8) of Kafiabad (2022) and can provide more accurate results.
Note that the evaluation of $\tilde{f}(\bm{Y}_{i},t_{n})$ requires an
interpolation since $\bm{Y}_{i}$ is not on the grid.

If we consider the particle position as a special case of $f$, (37) turns into

${\bm{\varXi}^{-1}}(\bm{X}_{i},t_{n+1})=\frac{(t_{n}-\tau)\bm{\varXi}^{-1}(\bm{Y}_{i},t_{n})+(t_{n+1}-t_{n})(\bm{X}_{i}+\bm{Y}_{i})/2}{t_{n+1}-\tau},$
(38)

which parallels Eq. (2.9) in Kafiabad (2022). The complete algorithm iterates
(36)–(38) from $t_{0}=\tau$ to $t_{N+1}=\tau+T$ to obtain ${\tilde{f}\mkern
1.5mu}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)=\tilde{f}(\bm{x},\tau+T)$
and $\bm{\varXi}^{-1}(\bm{x},t+T)$. In the final step,
$\bm{\varXi}^{-1}(\bm{x},t+T)$ is employed to apply the remapping (12) to
calculate ${\mathrel{\mathop{f}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{f}}}^{\scriptscriptstyle\mathrm{L}}(\bm{x},\tau)$.

##  Appendix C Lagrangian mean of tensors

We show how Eq. (17) for the Lagrangian mean of scalar functions generalises
readily to tensors when phrased in the language of push-forwards, pull-backs
and Lie derivatives (e.g. Frankel, 2004). The definition (6b) of the partial
Lagrangian mean generalises as

$\bar{\bm{\varphi}}^{{\scriptscriptstyle\mathrm{L}}*}_{\tau}{\mathrel{\mathop{\tau}\limits^{\vbox
to0.0pt{\kern-1.0pt\hbox{$\scriptstyle-\phantom{\mkern-5.5mu\rhook}$}\vss}}}}{\vphantom{\bar{\tau}}}^{\scriptscriptstyle\mathrm{L}}_{\tau}\coloneq\frac{1}{T}\int_{\tau}^{\tau+T}\bm{\varphi}_{s}^{*}\tau_{s}\,\mathrm{d}s.$
(39)

Here, we make the time variable appear explicitly as a subscript and we use
$\bar{\bm{\varphi}}^{{\scriptscriptstyle\mathrm{L}}*}_{t}$ and
$\bm{\varphi}_{t}^{*}$ to denote the pull-backs by, respectively, the
Lagrangian mean map and flow map. Pull-backs are interpreted differently
depending on the nature of $\tau_{t}$: in particular, for a scalar,
$(\bm{\varphi}_{t}^{*}f)(\bm{x})=f(\bm{\varphi}_{t}(\bm{x}))$, while for a
1-form,
$(\bm{\varphi}_{t}^{*}(\bm{u}\bcdot\mathrm{d}\bm{x}))(\bm{x})=u_{j}(\bm{\varphi}_{t}(\bm{x}))\partial_{i}\varphi_{tj}(\bm{x})\,\mathrm{d}x_{i}$.
The definition (39) is a natural, geometrically intrinsic definition of the
Lagrangian mean because it ensures that the tensors $\tau_{s}$ at different
times $s$ are transported to the same position in label space before the
averaging is carried out (see Gilbert & Vanneste, 2018). The partial mean of
$\tau_{t}$ associated with (39) reads

$\bar{\bm{\varphi}}_{t}^{*}\bar{\tau}_{t}\coloneq\frac{1}{t-\tau}\int_{\tau}^{t}\bm{\varphi}_{s}^{*}\tau_{s}\,\mathrm{d}s$
(40)

and is understood to depend parametrically on $\tau$.

Differentiating (40) with respect to $t$ and using that
$\partial_{t}\bar{\bm{\varphi}}_{t}^{*}\bar{\tau}_{t}=\bar{\bm{\varphi}}_{t}^{*}\left(\partial_{t}+\mathcal{L}_{\bar{\bm{u}}}\right)\bar{\tau}_{t}$,
where $\mathcal{L}_{\bar{\bm{u}}}$ is the Lie derivative, gives (Gilbert &
Vanneste, 2018)

$\bar{\bm{\varphi}}_{t}^{*}\left(\partial_{t}+\mathcal{L}_{\bar{\bm{u}}}\right)\bar{\tau}_{t}=\frac{-\bar{\bm{\varphi}}_{t}^{*}\bar{\tau}_{t}+\bm{\varphi}_{t}^{*}\bar{\tau}_{t}}{t-\tau}.$
(41)

Pushing forward by $\bar{\bm{\varphi}}_{t*}$ reduces this to

$\left(\partial_{t}+\mathcal{L}_{\bar{\bm{u}}}\right)\bar{\tau}_{t}=\frac{\bm{\varXi}_{t}^{*}\tau_{t}-\tau_{t}}{t-\tau}$
(42)

on using that (10), here in the form
$\bm{\varXi}_{t}=\bm{\varphi}_{t}\circ\bar{\bm{\varphi}}_{t}^{-1}$, implies
that $\bm{\varXi}_{t}^{*}=\bar{\bm{\varphi}}_{t*}\bm{\varphi}^{*}_{t}$. Eq.
(42) generalises (17) to tensors. Its coordinate form depends on the type of
tensor because of the different form taken by $\mathcal{L}_{\bar{\bm{u}}}$ and
$\bm{\varXi}_{t}^{*}$. Note that Eqs. (13) and (19) for $\bm{\varXi}^{-1}$ and
$\bm{\varXi}$ do not have a tensorial nature: because of the Cartesian
definition of the Lagrangian mean map (3), $\bm{\varXi}$ is simply a triple of
scalar functions rather than a vector. See Gilbert & Vanneste (2018) for
geometrically intrinsic definitions of the Lagrangian mean map alternative to
(3).

[Funding]HAK and JV were supported by the UK Engineering & Physical Sciences
Research Council (grant EP/W007436/1). JV was also 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 and scripts used for the simulations of
§4.2 are available at https://github.com/turbulencelover/ComputingLagMeans

[Author ORCID]H. A. Kafiabad, https://orcid.org/ 0000-0002-8791-9217; J.
Vanneste, https://orcid.org/0000-0002-0319-589X

## References

  * Andrews & McIntyre (1978) Andrews, D. G. & McIntyre, M. E. 1978 An exact theory of nonlinear waves on a Lagrangian-mean flow. J. Fluid Mech. 89, 609–646.
  * Bachman et al. (2020) Bachman, S. D., Shakespeare, C. J., Kleypas, J., Castruccio, F. S. & Curchitser, E. 2020 Particle-based Lagrangian filtering for locating wave-generated thermal refugia for coral reefs. J. Geophys. Res. Oceans 125, e2020JC016106.
  * Bretherton (1971) Bretherton, F. P. 1971 The general linearized theory of wave propagation. In Mathematical Problems in the Geophysical Sciences, Lect. Appl. Math., vol. 13, pp. 61–102. Am. Math. Soc.
  * Bühler (2014) Bühler, O. 2014 Waves and mean flows, 2nd edn. Cambridge University Press.
  * Bühler & McIntyre (1998) Bühler, O. & McIntyre, M. E. 1998 On non-dissipative wave–mean interactions in the atmosphere or oceans. J. Fluid Mech. 354, 301–343.
  * Frankel (2004) Frankel, T. 2004 The geometry of physics, 2nd edn. Cambridge University Press.
  * Gilbert & Vanneste (2018) Gilbert, A. D. & Vanneste, J. 2018 Geometric generalised Lagrangian-mean theories. J. Fluid Mech. 839, 95–134.
  * Gilbert & Vanneste (2021) Gilbert, A. D. & Vanneste, J. 2021 A geometric look at MHD and the Braginsky dynamo. Geophys. Astrophys. Fluid Dynam. 115, 436–471.
  * Grimshaw (1975) Grimshaw, R. 1975 Nonlinear internal gravity waves in a rotating fluid. J. Fluid Mech. 71 (3), 497–512.
  * Holm (2002) Holm, D. D. 2002 Lagrangian averages, averaged Lagrangians, and the mean effects of fluctuations in fluid dynamics. Chaos 12, 518–530.
  * Jones et al. (2022) Jones, C. S., Xiao, Q., Abernathey, R. & Smith, K. S. 2022 Separating balanced and unbalanced flow at the surface of the agulhas region using lagrangian filtering. EarthArXiv .
  * Kafiabad (2022) Kafiabad, H. A. 2022 Grid-based calculation of the Lagrangian mean. J. Fluid Mech. 940, A21.
  * Kafiabad et al. (2021) Kafiabad, H. A., Vanneste, J. & Young, W. R. 2021 Wave-averaged balance: a simple example. J. Fluid Mech. 911, R1.
  * Moore (1970) Moore, D. 1970 The mass transport velocity induced by free oscillations at a single frequency. Geophys. Fluid Dynam. 1, 237–247.
  * Nagai et al. (2015) Nagai, T., Tandon, A., Kunze, E. & Mahadevan, A. 2015 Spontaneous generation of near-inertial waves by the Kuroshio front. J. Phys. Oceanogr. 45 (9), 2381–2406.
  * Shakespeare et al. (2021) Shakespeare, C. J., Gibson, A. H., Hogg, A. McC., Bachman, S. D., Keating, S. R. & Velzeboer, N. 2021 A new open source implementation of Lagrangian filtering: A method to identify internal waves in high-resolution simulations. JAMES 13, e2021MS002616.
  * Shakespeare & Hogg (2017) Shakespeare, C. J. & Hogg, A. McC. 2017 Spontaneous surface generation and interior amplification of internal waves in a regional-scale ocean model. J. Phys. Oceanogr. 47, 811–826.
  * Shakespeare & Hogg (2018) Shakespeare, Callum J. & Hogg, Andrew McC. 2018 The life cycle of spontaneously generated internal waves. J. Phys. Oceanogr. 48, 343–359.
  * Shakespeare & Hogg (2019) Shakespeare, Callum J. & Hogg, Andrew McC. 2019 On the momentum flux of internal tides. J. Phys. Oceanogr. 49, 993–1013.
  * Soward (1972) Soward, A. M. 1972 A kinematic theory of large magnetic Reynolds number dynamos. Phil. Trans. R. Soc. London A272, 431–462.
  * Soward & Roberts (2010) Soward, A. M. & Roberts, P. H. 2010 The hybrid Euler–Lagrange procedure using an extension of Moffatt’s method. J. Fluid Mech. 661, 45–72.
  * Vallis (2017) Vallis, G. K. 2017 Atmospheric and oceanic fluid dynamics, 2nd edn. Cambridge University Press.
  * Vanneste (2013) Vanneste, J. 2013 Balance and spontaneous wave generation in geophysical flows. Annu. Rev. Fluid Mech. 45, 147–172.
  * Vanneste & Young (2022) Vanneste, J. & Young, W. R. 2022 Stokes drift and its discontents. Phil. Trans. R. Soc. A 380, 20210032.
  * Wagner & Young (2015) Wagner, G. L. & Young, W. R. 2015 Available potential vorticity and wave-averaged quasi-geostrophic flow. J. Fluid Mech. 785, 401–424.


# NASTAR: Noise Adaptive Speech Enhancement with Target-Conditional Resampling

###### Abstract

For deep learning-based speech enhancement (SE) systems, the training-test
acoustic mismatch can cause notable performance degradation. To address the
mismatch issue, numerous noise adaptation strategies have been derived. In
this paper, we propose a novel method, called noise adaptive speech
enhancement with target-conditional resampling (NASTAR), which reduces
mismatches with only one sample (one-shot) of noisy speech in the target
environment. NASTAR uses a feedback mechanism to simulate adaptive training
data via a noise extractor and a retrieval model. The noise extractor
estimates the target noise from the noisy speech, called pseudo-noise. The
noise retrieval model retrieves relevant noise samples from a pool of noise
signals according to the noisy speech, called relevant-cohort. The pseudo-
noise and the relevant-cohort set are jointly sampled and mixed with the
source speech corpus to prepare simulated training data for noise adaptation.
Experimental results show that NASTAR can effectively use one noisy speech
sample to adapt an SE model to a target condition. Moreover, both the noise
extractor and the noise retrieval model contribute to model adaptation. To our
best knowledge, NASTAR is the first work to perform one-shot noise adaptation
through noise extraction and retrieval.

Index Terms: speech enhancement, noise adaptation, contrastive learning,
source separation, acoustic retrieval

(a)

(b)

Figure 1: (a) Illustration of the subset relation between the accessible
source noise data set ($\mathcal{S}$), the target-similar noise set
$\tilde{\mathcal{T}}$, and the true target noise set ($\mathcal{T}$); (b) The
flowchart of our NASTAR system. Given a query signal, the noise extractor
directly estimates the pseudo-noise, and the retrieval model retrieves a
relevant-cohort set.

##  1 Introduction

Speech enhancement (SE) is commonly used as a pre-processor in speech-related
applications, such as hearing aids [1, 2], automatic speech recognition [3, 4,
5], and speech emotion recognition [6]. Recently, various deep learning (DL)
models have been used to formulate a regression function for SE [7, 8, 9, 10,
11, 12, 13]. Practically, we intend to collect a broad spectrum of noise types
to train a unified DL-based SE model that performs well in diverse conditions.
However, in real-world scenarios, the unified SE model inevitably encounters
unseen noise types. If the training data do not entirely align with the
testing distribution, a training-test mismatch occurs, limiting the achievable
enhancement performance.

Thus far, numerous attempts have been made to reduce the training-test
mismatch issue of SE. To prevent mismatch under specific noise types, noise-
aware strategies [14, 15, 16, 17] were proposed to specify noise information
into conditional embeddings to guide SE models to achieve better enhancement
results. Meanwhile, ensemble learning strategies [18, 19, 20] were used in SE
systems by preparing multiple component models, each of which is trained with
a subset of the training data; during testing, the outputs of these component
models are combined with a fusing/gating mechanism to dynamically handle noisy
conditions. Despite confirmed effectiveness in many tasks, noise-aware and
ensemble learning SE systems may yield sub-optimal performance when the
training data are not directly sampled for the target condition. To handle the
above issue, noise adaptation approaches have been proposed to handle
mismatches by fine-tuning SE models to properly match target conditions. For
supervised adaptation methods, the pre-trained SE system is modified with the
available noisy-clean paired signals [21, 22, 23, 24]. For unsupervised
adaptation methods, it is assumed that only noisy inputs are accessible.
Several studies [25, 26, 27] applied domain adversarial training (DAT) [28]
techniques to convert input speech signals to noise-invariant features for
processing by the decoder of the SE model. In addition, to avoid the
catastrophic forgetting problem, Lee et al. [29] proposed the SERIL system
based on a regularization-based incremental learning strategy. Meanwhile,
multi-task learning approaches [7, 30] employ auxiliary tasks to provide
positive effects for training adaptive SE models. All of these existing
strategies require batches of data for adaptation. However, we usually only
get a very small amount of data from the target environment.

**footnotetext: These authors contributed equally to this
work.$\dagger$$\dagger$footnotetext: The codes are availabel at
https://github.com/ChangLee0903/NASTAR

In this paper, we propose a novel method called noise adaptive speech
enhancement with target-conditional resampling (NASTAR†) to effectively handle
the mismatch problem with the least amount of target data. Instead of using
batches of extra data, our work is the first to reduce the mismatch problem by
filtering out out-of-target data. Given one segment of noisy speech in the
target environment, the NASTAR uses the sample to simulate the adaptation
training data via a noise extractor and a retrieval model. The noise extractor
estimates the target noise from the noisy speech, termed pseudo-noise. The
noise retrieval model retrieves a relevant noise set, termed relevant-cohort,
from an existing pool of noise signals based on the given noisy speech sample.
The pseudo-noise and the relevant-cohort set are used in a combined sampling
scheme and then mixed with the source speech corpus to simulate training data
for noise adaptation. Compared to existing approaches, the main advantages of
NASTAR include:

1) Data availability: NASTAR requires no additional target data beyond a given
“one-shot” from the target environment. In noise adaptive SE research, this is
the first work to re-utilize existing datasets instead of collecting
additional data.

2) Promising performance: Experiment results confirm that NASTAR yields
consistent and significant improvements (validated by the dependent t-test
shown in Sec. 4.2) over a pre-trained model with few speech samples.

3) High training stability: NASTAR performs adaptation in a simple supervised
manner, and the objective is committed to be aggressive towards handling the
target noise type. We note that NASTAR has a more stable training process than
the adaptive SE systems trained with multi-task objectives.

##  2 Background and Motivation

Based on the nature of noise interference, given a clean utterance $s(t)$ and
a noise signal $n(t)$, we can simulate the noisy signal by $s(t)+n(t)$ and
thus prepare noisy-clean paired data of the target condition. In general, we
assume that the recorded data from the target environment contain a mixture of
speech and noise components. Thus, we can use a noise extractor to estimate
the noise signal, termed pseudo-noise, from the noisy speech to simulate
target noisy speech for model adaptation. Nevertheless, the pseudo-noise may
not be sufficient to cover the environment characteristics if the target noise
is highly non-stationary, as there is only one segment. Therefore, we need to
leverage more target-similar noise signals for data simulation. Although the
prepared training data may cover a wide range of noise types, real-world SE
systems may encounter few specific noise types not fully involved in all the
training data, which inevitably leads to mismatches between training and test
data, thereby limiting the SE performance.

To address the mismatch issue, this study investigates retrieving similar
samples from the existing source noise dataset to enrich adaptation data. Our
main idea focuses on filtering out out-of-target data and finding “what to
train” for the specific target environment to reduce the mismatch. As shown in
Fig. 1(a), the proposed NASTAR system intends to collect the target-similar
noise set $\tilde{\mathcal{T}}$ from the source noise dataset $\mathcal{S}$
and the pseudo-noise signal to sufficiently overlap with the true target noise
set $\mathcal{T}$. Since our resampling method selects noise signals from
$\tilde{\mathcal{T}}$ instead of $\mathcal{S}$, this scheme is expected to
reduce the mismatch between training and target noisy environments. In this
way, the combination of noise extractor and noise retrieval model can provide
better noise adaptation results in a parallel training manner. Meanwhile, the
NASTAR system only needs one target noisy speech.

In our experimental setup, we firstly constructed a custom noise set, which
included $5$ different noise signals as different types. Each type of the
custom noise set corresponds to one target noise signal. We cut each target
noise signal in half. The first half was used to contaminate a randomly
selected test utterance as the accessible noisy speech (at SNR 0dB); the
second half was used to contaminate other test utterances.

##  3 The NASTAR System

###  3.1 System overview

Fig. 1(b) shows the overall flow of our proposed NASTAR system. First, we
assume that a segment of target noisy speech (one-shot) is available and
serves as the query signal. Next, the noise extractor estimates the noise from
the query signal as pseudo-noise. Meanwhile, the noise retrieval model
retrieves a relevant-cohort set from the source noise dataset based on the
query signal. In this work, the relevant-cohort set consists of 250 noise
signals closest to the query signal in terms of cosine similarity. The
collection of the pseudo-noise and the relevant-cohort set forms a noise pool
that potentially covers the noise characteristics of the target condition. By
using a collaborative sampling mechanism, the pseudo-noise signal has the
opportunity of “1 - $\alpha$”, and each noise signal in the relevant-cohort
set has the opportunity of “$\alpha$/250” to be selected to corrupt the clean
speech. Finally, we adapt the pre-trained SE model to the target condition by
using the sampled noise signals. $\alpha$ was set to 0.9.

###  3.2 Noise extractor

To construct a noise extractor, we used the DEMUCS model, which has been shown
to yield state-of-the-art results on the music separation and SE tasks [31,
32]. Referred to [32], we follow the same settings of the architecture and
objective function to estimate the noise signal as pseudo-noise from the query
signal. The DEMUCS model consists of an encoder-decoder architecture with
skip-connections and adopts a multi-resolution STFT objective function
capturing the information of different time-frequency resolutions more
effectively. The objective function is defined as:

$\frac{1}{T}[||y-\hat{y}||_{1}+\sum_{i=1}^{M}L_{stft}^{(i)}(y,\hat{y})],$ (1)

where $y$ and $\hat{y}$ are the target noise signal and the predicted noise
signal, respectively; $M$ is the number of STFT losses; $L_{stft}$ is the
addition of the spectral convergence loss $L_{sc}(y,\hat{y})$ and the
magnitude loss $L_{mag}(y,\hat{y})$. $L_{sc}(y,\hat{y})$ and
$L_{mag}(y,\hat{y})$ are respectively calculated by Eq. 2 and Eq. 3:

$L_{sc}(y,\hat{y})=\frac{\||STFT(y)|-|STFT(\hat{y})|\|_{F}}{\||STFT(y)|\|_{F}},$
(2) $L_{mag}(y,\hat{y})=\frac{1}{T}\|\log|STFT(y)|-\log|STFT(\hat{y})|\|_{1}.$
(3)

Figure 2: Formation of positive pairs and negative pairs for training the
noise retrieval model.

###  3.3 Noise retrieval model

Given a target noisy speech as the query signal, a noise retrieval model aims
to determine whether a compared noise signal is relevant to the query signal.
The key lies in a good representation of the query and compared signals. In
the representation space, the relevant signals are closer to the query than
the irrelevant signals. In addition, noise signals are diverse, complex, and
thus it is challenging to define a noise type as a specific class.

Some studies [33, 34] have explored self-supervised contrastive learning for
speaker verification as an unsupervised learning framework. In [34], the
commonly used contrastive learning frameworks SimCLR [35], MoCo [36], and
ProtoNCE [37] are compared for extracting speaker representations. To learn
discriminative noise representations, we refer to these frameworks and design
a pretext task to make paired data to train our noise retrieval model, as
shown in Fig. 2. Our pretext task used to produce paired segments from noise
signals include truncation, duplication, and shifting operations. Each sampled
noise signal is randomly processed into segments, and each segment ranges from
24,000 to 80,000 sample points. Since different segments in the same noise
signal have similar characteristics, these segments can be considered to be
relevant to each other. As shown in Fig. 2, given a segment (the query segment
$x_{q}$) excerpted from a noise signal, we can take another segment (the
positive segment $x_{k^{+}}$ to be compared) from the same noise signal to
form a positive pair. Meanwhile, we assume that any two noise signals come
from different noise conditions. Therefore, given a segment (the query segment
$x_{q}$), we can randomly select a segment from another noise signal (the
negative segment $x_{k^{-}}$ to be compared) to form a negative pair. In
addition, each segment being compared has a 50% chance of being mixed with a
speech signal at a random SNR level (from -8dB to 8dB in 2dB steps). With
these noisy speech segments involving both noise and speech components, our
noise retrieval model can learn to filter out non-noise (speech) components.
Therefore, at inference time, we can directly use “noisy speech” as the query
signal to find those acoustically similar noise samples to construct the
relevant-cohort set.

In training, we combine the negative pair usage mechanisms of the SimCLR [35]
and MoCo [36] frameworks. In the early stage, in order to stabilize the
training, we adopt the SimCLR method to generate negative pairs from the
segments in a batch. After 5,000 steps, we follow MoCo to use the queue to
access the past embeddings, and generate negative pairs from the past
embeddings and the current embeddings. For a mini-batch, there are two
embedding sets $\mathcal{B}_{q}$ and $\mathcal{B}_{k}$ generated by the query
encoder network $f_{q}$ and the key encoder network $f_{k}$, respectively, and
a queue $\mathcal{N}$ that stores embeddings in previous steps. The loss
function is defined as:

${l_{q,k^{+}}=-\log{\frac{\exp(\phi(q,k^{+})/\tau)}{\sum\limits_{k^{-}\in\mathcal{B}_{q}\cup\mathcal{B}_{k}\cup\mathcal{N}}\mathbbm{1}_{[k^{-}\neq
k^{+}]}\exp(\phi(q,k^{-})/\tau)}},}$ (4)

where $\phi$ is the similarity function, and $\tau$ denotes the temperature
parameter. The query embedding $q$ is obtained by $q=f_{q}(x_{q})$, and
$x_{q}$ is the query segment. The positive key embedding $k^{+}$ is obtained
by $k^{+}=f_{k}(x_{k}^{+})$, and $x_{k}^{+}$ is a positive key segment. For
the negative key embedding $k^{-}$ coming from $\mathcal{B}_{q}$ or
$\mathcal{B}_{k}$, $k^{-}$ is obtained by $k^{-}=f_{q}(x_{k}^{-})$ or
$k^{-}=f_{k}(x_{k}^{-})$, and $x_{k}^{-}$ is a negative key segment. The
negative key embedding $k^{-}$ can also be sampled from the previous
embeddings $\mathcal{N}$. Given an embedding pair $(q,k)\in\mathbb{R}^{d}$, we
use cosine similarity in the similarity function $\phi$.

Our noise retrieval model consists of three bidirectional LSTM layers and one
cascaded MLP projection layer. Denoting the parameters of $f_{k}$ as
$\theta_{k}$ and those of $f_{q}$ as $\theta_{q}$, we update $\theta_{k}$ by
$\theta_{k}=\mu\theta_{k}+(1-\mu)\theta_{q}$, where $\mu$ is the momentum
coefficient. To stabilize training in the early stage, those negative keys
stored in $\mathcal{N}$ are not sampled until 5,000 steps.

Table 1: Average $\text{PESQ}_{\text{nb}}$, STOI, and SI-SDR scores of different models. The red and orange numbers denote the first and second place results for each condition, respectively. | ACVacuum | Babble | CafeRestaurant | Car | MetroSubway
---|---|---|---|---|---
method | STOI |  PESQ${}_{\text{nb}}$ | SI-SDR | STOI |  PESQ${}_{\text{nb}}$ | SI-SDR | STOI |  PESQ${}_{\text{nb}}$ | SI-SDR | STOI |  PESQ${}_{\text{nb}}$ | SI-SDR | STOI |  PESQ${}_{\text{nb}}$ | SI-SDR
NOISY | 0.8407 | 1.9671 | 10.8404 | 0.8469 | 2.2791 | 12.3311 | 0.8955 | 2.7286 | 10.8916 | 0.8706 | 2.4135 | 11.9138 | 0.9359 | 3.3332 | 21.2817
PTN | 0.8815 | 2.8496 | 17.6623 | 0.8862 | 2.8116 | 16.4004 | 0.9199 | 3.3711 | 17.5752 | 0.925 | 3.4061 | 16.8882 | 0.9488 | 3.8564 | 21.2581
DAT${}_{\text{one}}$ | 0.8823 | 2.8694 | 17.8753 | 0.8897 | 2.8017 | 16.9236 | 0.9215 | 3.3794 | 17.4182 | 0.9276 | 3.4170 | 17.2144 | 0.9504 | 3.8169 | 22.3700
DAT${}_{\text{full}}$ | 0.8712 | 2.5442 | 15.5812 | 0.884 | 2.7615 | 16.1574 | 0.9094 | 3.2511 | 16.7324 | 0.9089 | 3.0483 | 15.9469 | 0.9359 | 3.7047 | 20.5600
NASTAR | 0.8929 | 2.9482 | 18.2684 | 0.8916 | 2.8655 | 17.2072 | 0.9244 | 3.428 | 18.6476 | 0.930 | 3.4602 | 18.0286 | 0.9524 | 3.9121 | 22.3703
EXTR | 0.8914 | 2.9047 | 18.6487 | 0.8771 | 2.6580 | 15.9827 | 0.9144 | 3.2114 | 14.9499 | 0.9113 | 3.1036 | 17.3962 | 0.9439 | 3.7233 | 21.5684
GT | 0.8929 | 2.9062 | 18.7776 | 0.8779 | 2.6576 | 16.1166 | 0.9167 | 3.2485 | 16.3388 | 0.9106 | 3.0888 | 17.2334 | 0.9457 | 3.7811 | 23.3171
ALL | 0.8905 | 2.9330 | 18.2512 | 0.8903 | 2.8260 | 16.8602 | 0.9225 | 3.4186 | 18.4044 | 0.9287 | 3.4225 | 17.4977 | 0.9504 | 3.8201 | 22.0051
RETV | 0.8883 | 2.9469 | 18.1033 | 0.8913 | 2.8672 | 17.1955 | 0.9244 | 3.438 | 18.3673 | 0.9282 | 3.4949 | 18.2736 | 0.9524 | 3.8879 | 22.7825
OPT | 0.9017 | 3.0561 | 18.9691 | 0.9111 | 3.0662 | 18.5652 | 0.9348 | 3.5684 | 19.4992 | 0.9364 | 3.5641 | 19.4694 | 0.9568 | 3.9346 | 26.5219

##  4 Experiments

###  4.1 Experimental setup and implementation details

The datasets used in the experiments include: (1) the DNS-Challenge [38] noise
dataset consisting of 65,303 background and foreground noise samples; (2) the
Librispeech-360 [39] corpus consisting of 104,014 utterances; (3) the Voice
Bank Corpus (VBC) [39] consisting of 11,572 utterances of 28 speakers; and (4)
the custom noise set consisting of five noise signals of five common noise
types, namely AC/Vacuum, Babble, CafeRestaurant, Car, and MetroSubway,
selected from FreeSound [40].

The NASTAR system consists of four modules: the noise extractor, the noise
retrieval model, the pre-trained SE model, and the adapted SE model. To train
the noise extractor and the pre-trained SE model, the clean speech utterances
from Librispeech-360 were contaminated by randomly sampled noise signals from
DNS-Challenge at five SNR levels (from 0dB to 12dB in 3dB steps). We prepared
noisy-noise and noisy-speech pairs to train the noise extractor and the pre-
trained SE model, respectively. The models were trained by the Adam optimizer
with $\beta_{1}$ = 0.9 and $\beta_{2}$ = 0.999, a learning rate of 0.0002, and
a batch size of 8 with 500,000 steps. The noise extractor and the pre-trained
SE model shared the same architecture and hyper-parameter settings. To train
the noise retrieval model, the training data were prepared from DNS-Challenge
and Librispeech-360. The process of preparing positive and negative pairs is
presented in Sec. 3.3 and Fig. 2. The noise retrieval model was trained by
contrasting positive and negative pairs. We used the Adam optimizer with
$\beta_{1}$ = 0.9 and $\beta_{2}$ = 0.999, a learning rate of 0.00025, and a
batch size of 256 to train the model for 100,000 steps. The hyper-parameters
were set as follows, the temperature coefficient $\tau$ = 0.1, the momentum
coefficient $\mu$ = 0.9, and the queue size $|\mathcal{N}|$ = 32,768. To train
the adapted SE model, the pseudo-noise and the relevant-cohort set were
derived based on the given query signal. The pre-trained SE model served as
the initial checkpoint. Each clean utterance from Librispeech-360 was mixed
with a target-similar noise signal sampled from the pseudo-noise and relevant-
cohort at 7 SNR levels (from -4dB to 8dB in 2dB steps). The SE model was
updated by the Adam optimizer with $\beta_{1}$ = 0.9 and $\beta_{2}$ = 0.999,
a learning rate of 0.0001, and a batch size of 8 with 20,000 steps.

For each noise signal in the custom noise set, we used its first half and a
test speech from the test-clean set of Librispeech to form the 0dB query noisy
speech and accordingly trained the adaptive SE model. That is, there are five
adaptation models, each for one noise type. For evaluation, we used the clean
speech of the VBC corpus as our test set. For each noise condition, each clean
utterance was separately mixed with the second half of the corresponding
target noise signal at four SNR levels (from -5dB to 10dB in 5dB steps) as the
test noisy speech. Therefore, there are 11,572$\times$4 test noisy utterances
for each noise condition. Three standardized evaluation metrics were used to
measure the SE performance: narrow band perceptual evaluation of speech
quality ($\text{PESQ}_{\text{nb}}$) [41], short-time objective intelligibility
measure (STOI) [42], and scale-invariant signal-to-distortion ratio (SI-SDR)
[43]. $\text{PESQ}_{\text{nb}}$ was designed to evaluate the quality of
processed speech, ranging from -0.5 to 4.5. STOI was designed to compute the
speech intelligibility, ranging from 0 to 1. SI-SDR was designed to measure
the energy ratio between speech and non-speech components. For all three
metrics, higher scores indicate better performance.

###  4.2 Results

Table 1 shows the STOI, $\text{PESQ}_{\text{nb}}$, and SI-SDR scores of NASTAR
and other related SE methods. PTN represents the original pre-trained SE
model. DAT${}_{\text{one}}$ represents the adaptive SE model trained with DAT
using one sample of the test set. DAT${}_{\text{full}}$ represents the
adaptive SE model trained with DAT using the full test set. To investigate the
effect of the noise extractor and the noise retrieval model in our NASTAR
model, we prepared the paired data for adaptation in different ways: (1) only
based on the estimated pseudo-noise (termed EXTR, $\alpha$ = 0); (2) only
based on the ground truth of pseudo-noise (termed GT, $\alpha$ = 0); (3) based
on the estimated pseudo-noise and randomly sampled noise signals from DNS-
Challenge (termed ALL, $\alpha$ = 0.9); and (4) only based on the 250 noise
signals in relevant-cohort (termed RETV, $\alpha$ = 1). In addition, we used
the second half of the target noise signal to adapt the pre-trained SE model;
that is, the adaptation and test noise signal are the same, and $\alpha$ is 0.
The resulting model is called OPT. The oracle performance by OPT represents
the upper-bound of NASTAR’s performance. Furthermore, our test set can be
divided into 20 groups based on 5 noise types and 4 SNR levels. We performed a
dependent t-Test on the paired SE results based on these 20 groups. The
p-values for NASTAR to outperform PTN in STOI, PESQ${}_{\text{nb}}$, and SI-
SDR scores are $3.31\text{\times}{10}^{-8}\text{\,}\mathrm{,}$
$5.81\text{\times}{10}^{-11}\text{\,}\mathrm{,}$ and
$3.35\text{\times}{10}^{-12}\text{\,}\mathrm{,}$ respectively. Based on a
standard threshold of 0.05, the improvement of NASTAR over PTN is confirmed to
be significant.

####  4.2.1 Comparison with other adaptive models

Next, we compare the performance of NASTAR and the adaptive SE model with
domain adversarial training (cf. DAT${}_{\text{one}}$ and
DAT${}_{\text{full}}$ in Table 1). From Table 1, we observe that
DAT${}_{\text{one}}$ and DAT${}_{\text{full}}$ are not always better than PTN,
and our NASTAR model consistently outperforms DAT${}_{\text{one}}$ and
DAT${}_{\text{full}}$ in all metrics under all noise conditions. The main
reason for the results is that the help of the auxiliary adversarial training
loss cannot effectively improve SE even with complete test data. In contrast,
our NASTAR model is trained to the original SE objective by adjusting the
sampling scheme of the existing training data.

####  4.2.2 Ablation study

Then, we investigate the effect of the noise extractor and the noise retrieval
model in NASTAR. From Table 1, we observe that EXTR performs mostly worse than
other variants of the NASTAR model. The results show that due to the
estimation error of the noise extractor, the low-quality pseudo-noise is not
sufficient to cover the environmental characteristics, and the complement of
relevant-cohort is needed. Comparing NASTAR with ALL, we confirm that adaptive
data sampled from relevant-cohort are more useful than randomly sampled data
from DNS-Challenge. This is because the noise samples in relevant-cohort are
closer to the target condition than the randomly sampled data from DNS-
Challenge, thus reducing the mismatch between training and test conditions. GT
performs mostly worse than NASTAR and RETV, from which we can see that the
sample diversity is not enough to cope with the variation of target noise. In
contrast, the noise samples in relevant-cohort used in NASTAR and RETV are
close to the target condition, thus enriching the adaptive data. Furthermore,
GT is not always better than EXTR. This result shows that the performance of
the noise extractor is acceptable. Overall, NASTAR performs the best, followed
by RETV. The results confirm the effectiveness of the noise extractor and
noise retrieval model in NASTAR. NASTAR is slightly worse than OPT. There is
still room for improvement.

Figure 3: The average relative improvement rate in STOI and
PESQ${}_{\text{nb}}$ scores for different SE methods.

####  4.2.3 Relative improvement rate

From Table 1, we note that the STOI improvement of the pre-trained SE model
(PTN) over noisy speech is 0.0408 (0.8815 - 0.8407) in AC/Vacuum, and the
improvement of NASTAR over PTN is 0.0114 (0.8929 - 0.8815) in AC/Vacuum. We
also note that NASTAR has a significant improvement over PTN in terms of
$\text{PESQ}_{\text{nb}}$ and SI-SDR scores. As for the other four noise
types, we observe the same results that NASTAR consistently outperforms PTN in
terms of $\text{PESQ}_{\text{nb}}$, STOI, and SI-SDR scores. To further
evaluate the improvement brought by adaptation, we define a relative
improvement rate as
$\frac{s-s_{\text{NOISY}}}{s_{\text{PTN}}-s_{\text{NOISY}}}$, where $s$ is the
metric score obtained by an adaptation method, and $s_{\text{NOISY}}$ is the
metric score of the unprocessed noisy speech, and $s_{\text{PTN}}$ is the
metric score obtained by the pre-trained SE model (PTN). The relative
improvement rate indicates how much further improvement the adaptation method
can provide compared to PTN. Fig. 3 shows the average relative improvement
rate in STOI and PESQ${}_{\text{nb}}$ scores for NASTAR and other SE systems
across five noise types. It is clear from the figure that NASTAR reached top
performance in both STOI and PESQ${}_{\text{nb}}$ compared to other SE
methods, confirming the effectiveness of NASTAR.

##  5 Concluding Remarks

Existing noise adaptation strategies for SE require batches of target noisy
speech signals; however, only one-shot can be obtained in most cases. In this
paper, we have proposed a novel NASTAR system to reduce the acoustic mismatch
with only one target noisy utterance. NASTAR uses an accessible sample to
simulate adaptation data via a noise extractor and a noise retrieval model.
The noise extractor estimates the target noise, termed pseudo-noise, in the
noisy speech sample. The noise retrieval model retrieves a relevant noise set,
termed relevant-cohort, from a noise signal pool according to the given noisy
speech sample. The pseudo-noise and relevant-cohort set are sampled and mixed
with the source speech corpus to prepare simulated data for SE model
adaptation. Experimental results show that NASTAR can effectively utilize a
single target sample and outperforms several existing methods under different
target conditions. To our best knowledge, NASTAR is the first work to apply
resampling to leverage existing source data, which is more precise, efficient,
and effective. Furthermore, our adaptation method can be combined with
existing noise-aware and ensemble strategies in various scenarios.

## References

  * [1] K. Tan, X. Zhang, and D. Wang, “Real-time speech enhancement using an efficient convolutional recurrent network for dual-microphone mobile phones in close-talk scenarios,” in _Proc. of ICASSP_ , 2019.
  * [2] I. Fedorov, M. Stamenovic, C. Jensen, L. Yang, A. Mandell, Y. Gan, M. Mattina, and P. N. Whatmough, “Tinylstms: Efficient neural speech enhancement for hearing aids,” in _Proc. of Interspeech_ , 2020.
  * [3] A. Nicolson and K. K. Paliwal, “Deep Xi as a front-end for robust automatic speech recognition,” in _Proc. of CSDE_ , 2020.
  * [4] Z. Meng, S. Watanabe, J. R. Hershey, and H. Erdogan, “Deep long short-term memory adaptive beamforming networks for multichannel robust speech recognition,” in _Proc. of ICASSP_ , 2017.
  * [5] F. Weninger, H. Erdogan, S. Watanabe, E. Vincent, J. Le Roux, J. R. Hershey, and B. Schuller, “Speech enhancement with lstm recurrent neural networks and its application to noise-robust asr,” in _Latent Variable Analysis and Signal Separation_ , E. Vincent, A. Yeredor, Z. Koldovský, and P. Tichavský, Eds. Springer International Publishing, 2015, pp. 91–99.
  * [6] A. Triantafyllopoulos, G. Keren, J. Wagner, I. Steiner, and B. W. Schuller, “Towards robust speech emotion recognition using deep residual networks for speech enhancement,” in _Proc. of Interspeech_ , 2019.
  * [7] S. Wang, W. Li, S. M. Siniscalchi, and C. Lee, “A cross-task transfer learning approach to adapting deep speech enhancement models to unseen background noise using paired senone classifiers,” in _Proc. of ICASSP_ , 2020.
  * [8] Y.-C. Lin, Y.-T. Hsu, S.-W. Fu, Y. Tsao, and T.-W. Kuo, “IA-NET: Acceleration and compression of speech enhancement using integer-adder deep neural network,” in _Proc. of Interspeech_ , 2019.
  * [9] X. Lu, Y. Tsao, S. Matsuda, and C. Hori, “Speech enhancement based on deep denoising autoencoder.” in _Proc. of Interspeech_ , 2013.
  * [10] Y. Xu, J. Du, L.-R. Dai, and C.-H. Lee, “A regression approach to speech enhancement based on deep neural networks,” _IEEE/ACM TASLP_ , vol. 23, no. 1, pp. 7–19, 2014.
  * [11] D. Liu, P. Smaragdis, and M. Kim, “Experiments on deep learning for speech denoising,” in _Proc. of Interspeech_ , 2014.
  * [12] H. Erdogan, J. R. Hershey, S. Watanabe, and J. Le Roux, “Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks,” in _Proc. of ICASSP_ , 2015.
  * [13] A. Kumar and D. Florencio, “Speech enhancement in multiple-noise conditions using deep neural networks,” in _Proc. of Interspeech_ , 2016.
  * [14] “Noise perturbation for supervised speech separation,” _Speech Communication_ , vol. 78, pp. 1–10, 2016.
  * [15] B. O. Odelowo and D. V. Anderson, “A study of training targets for deep neural network-based speech enhancement using noise prediction,” in _Proc. of ICASSP_ , 2018.
  * [16] H. Li and J. Yamagishi, “Noise Tokens: Learning Neural Noise Templates for Environment-Aware Speech Enhancement,” in _Proc. of Interspeech_ , 2020\.
  * [17] J. Lee, Y. Jung, M. Jung, and H. Kim, “Dynamic noise embedding: Noise aware training and adaptation for speech enhancement,” in _Proc. of APSIPA ASC_ , 2020.
  * [18] X.-L. Zhang and D. Wang, “A deep ensemble learning method for monaural speech separation,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 2016.
  * [19] S. E. Chazan, J. Goldberger, and S. Gannot, “Deep recurrent mixture of experts for speech enhancement,” in _Proc. of WASPAA_ , 2017.
  * [20] J. Le Roux, S. Watanabe, and J. R. Hershey, “Ensemble learning for speech enhancement,” in _2013 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics_ , 2013.
  * [21] Y. Xu, J. Du, L.-R. Dai, and C.-H. Lee, “Cross-language transfer learning for deep neural network based speech enhancement,” in _Proc. of ISCSLP_ , 2014\.
  * [22] D. Wang and T. F. Zheng, “Transfer learning for speech and language processing,” in _Proc. of APSIPA ASC_ , 2015.
  * [23] Y. Koizumi, K. Niwa, Y. Hioka, K. Kobayashi, and Y. Haneda, “Dnn-based source enhancement to increase objective sound quality assessment score,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 2018.
  * [24] ——, “Dnn-based source enhancement self-optimized by reinforcement learning using sound quality measurements,” in _Proc. of ICASSP_ , 2017.
  * [25] C.-F. Liao, Y. Tsao, H.-Y. Lee, and H.-M. Wang, “Noise adaptive speech enhancement using domain adversarial training,” in _Proc. of Interspeech_ , 2019.
  * [26] T. Higuchi, K. Kinoshita, M. Delcroix, and T. Nakatani, “Adversarial training for data-driven speech enhancement without parallel corpus,” in _Proc. of IEEE ASRU_ , 2017.
  * [27] H.-Y. Lin, H.-H. Tseng, X. Lu, and Y. Tsao, “Unsupervised noise adaptive speech enhancement by discriminator-constrained optimal transport,” in _Proc. of NeurIPS_ , 2021.
  * [28] Y. Ganin, E. Ustinova, H. Ajakan, P. Germain, H. Larochelle, F. Laviolette, M. Marchand, and V. Lempitsky, “Domain-adversarial training of neural networks,” _JMLR_ , vol. 17, no. 1, pp. 2096–2030, 2016.
  * [29] C.-C. Lee, Y.-C. Lin, H.-T. Lin, H.-M. Wang, and Y. Tsao, “SERIL: noise adaptive speech enhancement using regularization-based incremental learning,” in _Proc. of Interspeech_ , 2020.
  * [30] Y. Bando, K. Sekiguchi, and K. Yoshii, “Adaptive neural speech enhancement with a denoising variational autoencoder.” in _Proc. of Interspeech_ , 2020, pp. 2437–2441.
  * [31] A. Défossez, N. Usunier, L. Bottou, and F. Bach, “Music source separation in the waveform domain,” _arXiv preprint arXiv:1911.13254_ , 2019.
  * [32] A. Defossez, G. Synnaeve, and Y. Adi, “Real time speech enhancement in the waveform domain,” in _Proc. of Interspeech_ , 2020.
  * [33] M. Ravanelli and Y. Bengio, “Learning Speaker Representations with Mutual Information,” in _Proc. of Interspeech_ , 2019.
  * [34] W. Xia, C. Zhang, C. Weng, M. Yu, and D. Yu, “Self-supervised text-independent speaker verification using prototypical momentum contrastive learning,” in _Proc. of ICASSP_ , 2021.
  * [35] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” in _Proc. of ICML_ , 2020\.
  * [36] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” in _Proc. of CVPR_ , 2020.
  * [37] J. Li, P. Zhou, C. Xiong, and S. Hoi, “Prototypical contrastive learning of unsupervised representations,” in _Proc. of ICLR_ , 2021.
  * [38] C. K. Reddy, V. Gopal, R. Cutler, E. Beyrami, R. Cheng, H. Dubey, S. Matusevych, R. Aichner, A. Aazami, S. Braun _et al._ , “The interspeech 2020 deep noise suppression challenge: Datasets, subjective testing framework, and challenge results,” in _Proc. of Interspeech_ , 2020\.
  * [39] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “Librispeech: an asr corpus based on public domain audio books,” in _Proc. of ICASSP_ , 2015.
  * [40] F. Font, G. Roma, and X. Serra, “Freesound technical demo,” in _Proc. of the 21st ACM International Conference on Multimedia_ , 2013.
  * [41] A. W. Rix, J. G. Beerends, M. P. Hollier, and A. P. Hekstra, “Perceptual evaluation of speech quality (PESQ)-a new method for speech quality assessment of telephone networks and codecs,” in _Proc. of ICASSP_ , 2001.
  * [42] C. Taal, R. Hendriks, R. Heusdens, and J. Jensen, “A short-time objective intelligibility measure for time-frequency weighted noisy speech,” in _Proc. of ICASSP_ , 2010.
  * [43] J. Le Roux, S. Wisdom, H. Erdogan, and J. Hershey, “SDR-half-baked or well done?” in _Proc. of ICASSP_ , 2019.