Volume-preserving Lagrangian averaging using polar factorization

Introduction

Averaging is a key tool for the mathematical and numerical modelling of fluid flows with multiple time scales. Standard averaging is Eulerian, that is, carried out at fixed spatial location. An alternative is Lagrangian averaging, carried out at fixed particle label, that is, along particle trajectories. Lagrangian averaging has conceptual and practical advantages over Eulerian averaging. Conceptually, Lagrangian averaging leads to averaged dynamical equations that are often simpler and closer in structure to the original unaveraged equations than those obtained with Eulerian averaging [?, ?, ?, ?, 5, 6, 7, 8, 9]. Practically, Lagrangian averaging avoids the entanglement between fast and slow motion that results from their mutual advection. In particular, it eliminates the blurring of flow features that fast advection introduces in Eulerian averages. This enables a clean separation of the flow into slow and fast components each often associated with vortical and wave motion [10, 11, 12, 13, 14].

The generalised Lagrangian mean (GLM) theory of Andrews and McIntyre [?] provides a framework for the formulation of Lagrangian-averaged models. A key feature is the replacement of particle labels by Lagrangian mean positions as independent variables for the Lagrangian mean fields. The Lagrangian mean positions are usually defined via a straightforward component-wise average of the actual positions. While this definition has the advantage of simplicity, it introduces the so-called divergence effect [9, 15]: for an incompressible fluid, the Lagrangian mean velocity is divergent even though the actual velocity is divergence free. Equivalently, the map between mean and actual positions does not preserve volume.

To illustrate the divergence effect, consider a two-dimensional incompressible fluid in rigid-body rotation. The trajectories $\boldsymbol{x}(t) = (x(t),y(t))$ of particles are circles. It is intuitively clear that the component-wise average

where $T$ is the averaging time, defines a mean position $\bar{\boldsymbol{x}}(t)$ with

see Figure 1. Thus time averaging contracts circular trajectories. A simple computation in appendix A confirms this conclusion and extends it to a broad class of time averages.

Alternatives to the component-wise averaging of positions have been proposed to remedy the divergence effect and obtain volume-preserving versions of Lagrangian averaging [?, 16, 17, 18]. These have been introduced for perturbative analyses of wave–mean flow process. In this context, they can be regarded as small modifications of standard GLM that trade the simplicity afforded by volume preservation (a constant mean density) against other complications (the perturbation fields acquire non-zero means).

In this paper, we revisit volume-preserving Lagrangian averaging from a different perspective. The new perspective is brought about by the recent development of efficient, PDE-based methods for the numerical computation of Lagrangian averages from simulation data [13, 19, 20, 21]. These methods, used primarily as diagnostic tools, rely on the standard GLM definition of the mean position and therefore suffer from the divergence effect. Here we overcome this drawback by developing a numerical method for the computation of a volume-preserving Lagrangian average.

The method is based on the idea put forward by Gilbert and Vanneste [17] of replacing the mean flow map (sending particle labels to mean positions) of GLM by the closest volume-preserving map. Closest is here understood in a straightforward mean-square sense. With this definition, we can exploit a seminal result of optimal transport theory [22] and characterize the new mean map as the volume-preserving factor in the polar factorization of the GLM map.

We obtain a practical formulation, based on the velocity of the flow rather than the (usually exceedingly intricate) flow map, by considering the time derivative of the polar factorization. This leads to a system of equations relating the divergence-free Lagrangian mean velocity to the flow velocity (§II). This system is similar to that obtained for GLM in that it couples the Lagrangian mean velocity with the map between mean and actual positions, but it also involves an additional scalar field that appears in the polar factorization.

The construction of a volume-preserving mean map via polar factorization is not perturbative, that is, it does not require the common assumption that the flow consists of a large mean component plus small-amplitude perturbations. Under this assumption, we recover results of Ref. [17] and clarify the differences between the mean flows corresponding to GLM, to polar factorization, and to the alternative volume-preserving construction in Refs. [16, 18] (appendix C). We also show how a superficially appealing solution to the divergence effect – the projection of the GLM velocity onto a divergence-free velocity – fails because it leads to mean trajectories that drift away from the particle trajectories.

We demonstrate the capability of the proposed volume-preserving formulation for the computation of Lagrangian time averages from simulation data (§III). We choose the exponential filter as averaging operator. This is the simplest filter in a family of sum-of-exponential filters which make it possible to compute averages alongside the dynamical equations in a truly on-the-fly fashion, through the solution of additional evolution equations [21]. The volume-preserving implementation retains much of the simplicity of the original GLM implementation in Ref. [21], requiring only a modest additional computational effort.

We implement the volume-preserving exponential filter for simulations of a two-dimensional incompressible fluid (§III B) and of a rapidly rotating shallow water system (§III C). We use the incompressible case, in which there are no fast waves to filter out, to contrast the way in which GLM and volume-preserving averaging alter slow vertical motion. The rotating shallow-water flow superposes to a turbulent vortical flow a fast Poincaré wave which the Lagrangian averaging aims to filter out. Although the shallow-water model is compressible, we find benefits in imposing volume-preservation of the mean map in the form of better wave filtering. We relate this finding to the approximate geostrophic balance of Lagrangian mean flows in the rapidly rotating regime [e.g. [23]].

Applications of the central idea of this paper, namely the polar factorization of the GLM mean map, are not restricted to the exponential filter. We exemplify this with an implementation of the 2nd-order Butterworth filter, a more advanced sum-of-exponential filter (appendix D). We expect that the idea can also be adapted to more general filters used for online or offline computations of volume-preserving Lagrangian means.

Formulation

GLM

The motion of a fluid is determined by the flow map $\varphi(a,t)$ which gives the position at time $t$ of particles labeled by $a$ in a reference configuration. Lagrangian averaging starts by defining an analogous mean flow map $\bar{\varphi}(a,t)$. In GLM, it is simply the component-wise average of the flow map:

for $i=1,\,2,\,3$. Since we restrict our attention to Euclidean space, we can treat maps as vectors and rewrite (3) as

In what follows we consider two types of averages: time averages of the form

where the kernel satisfies

and ensemble averages of the form

where $\varpi$ is a ensemble parameter with probability density $k(\varpi)$ such that

Ensemble averages assume an ensemble of flow maps, $\varphi(a,t,\varpi)$; we drop the explicit dependence on $\varpi$ in what follows, but it should be kept in mind.

The Lagrangian mean velocity, which we denote by $\bar{u}$ (without a superscript L, to keep the notation light), is the velocity associated with the mean flow map:

With (4), it is also the (componentwise) Lagrangian mean of the velocity in the sense that

The flow map is decomposed as

where $\Xi$ is the perturbation map, which depends on $\varpi$ for ensemble averaging. With this decomposition and for ensemble average the definition (10) of the Lagrangian mean velocity can be rewritten as

by substituting $x$ for the mean position, $x=\bar{\varphi}(a,t)$. The same substitution in (4) gives

With $\Xi(x,t)=x+\xi(x,t)$, this becomes

Eqs. (12) and (14) are taken as defining equations for GLM by Andrews and McIntyre [?]. They do not hold for time averages, however, because of the time dependence of the substitution $x=\bar{\varphi}(a,t)$. As a result, (12) and (14) only hold approximately provided that there is a time-scale separation between mean and perturbation.

Volume preservation

Mass conservation implies that the density $\rho(x,t)$ of a fluid is related to a fixed density $\rho_0(a)$ in label space according to

where $x=\varphi(a,t)$. Expressing the volume element $\mathrm{d}x$ in terms of $\mathrm{d}a$ then gives

Here $\left|\partial\varphi/\partial a\right|$ with $\lvert\cdot\rvert$ denoting determinant is the Jacobian of the flow map. The continuity equation

follows by time differentiation and use of the Jacobi identity.

For a homogeneous incompressible fluid, the densities $\rho$ and $\rho_0$ are constants which can be taken to be 1: $\rho=\rho_0=1$. The flow map $\varphi(a,t)$ preserves volume and its Jacobian is unity. Because this constraint is nonlinear, its average is non-trivial and the mean map does generally not preserve volume. Hence,

Correspondingly, while the velocity of an incompressible fluid is divergence-free, the Lagrangian mean velocity is divergent,

This is the divergence effect in GLM. It necessitates the introduction of a density variable $\widetilde{\rho}(x,t)$ with

or, equivalently since $\rho^{-1}=\left|\partial\varphi/\partial a\right|=\left|\partial\Xi/\partial x\right|\left|\partial\bar{\varphi}/\partial a\right|=1$,

which satisfies the continuity equation

analogous to (17).

It is possible to modify the definition of the mean flow map, hence (12) and (14), to eliminate the divergence effect. We propose to replaces the “bare” GLM mean map (4) by the closest flow map that preserves volume, that is, to replace $\bar{\varphi}$ by

where the integration is over the reference configuration. Here and in what follows we use the subscript $\dagger$ to indicate fields associated with the volume-preserving construction. In (23) we measure closeness in terms of an $L_2$ norm; other choices are possible.

A definition equivalent to (23) for ensemble averaging is

as introduced by Gilbert and Vanneste [17]. It states that $\bar{\varphi}^{\dagger}$ is the Fréchet mean of the maps $\bar{\varphi}$ for the $L_2$ distance. The equivalence follows from the equalities

using (4). The equivalence does not apply to time averaging because $\psi(a,t)\cdot\bar{\varphi}(a,t)\ne\langle\psi(a,t)\cdot\varphi(a,t)\rangle$.

The minimization (23) is a standard problem in the theory of optimal transport [e.g. [24]]. It leads to Brenier's [22] polar factorization: the bare mean map can be written as the composition

of the gradient of a convex scalar function $\lambda$ with the desired volume-preserving flow map. The characterization (26) is central to the remainder of the paper. For completeness, we provide a formal derivation in appendix B.

Given $\bar{\varphi}$, it is in principle possible to determine $\bar{\varphi}^{\dagger}$ by solving either the minimization problem (23) or a Monge–Ampère equation for the Legendre transform of $\lambda$ [22]. For the purpose of computing Lagrangian means from simulation data, it is better to consider the time derivative of (26) to obtain volume-preserving versions of (12) and eq:14`.

Differentiating (26) with respect to $t$ and using (10) gives

The Lagrangian mean velocity $\bar{u}^{\dagger}$ associated with $\bar{\varphi}^{\dagger}$ satisfies

and is divergence free:

Using (28) we rewrite (27) as

For the ensemble average we can freely move $\bar{\varphi}^{\dagger}$ on the right-hand side within the average. Using the decomposition

then gives

This equation is conveniently replaced by its divergence and curl, that is,

On the other hand, differentiating (31) with respect to time leads to

Eqs. (29), (33), (34) and (35) are key results of this section. They constitute a closed system that determines $\bar{u}^{\dagger}$, $\lambda$ and $\Xi^{\dagger}$ from $u$. Specifically, (33) and (35) govern the joint evolution of $\lambda$ and $\Xi^{\dagger}$ while $\bar{u}^{\dagger}$ is diagnosed from (29) and (34). This system replaces the simpler GLM prescription that couples (35) (without the dagger subscripts) and the relation $\bar{u}=\langle u\circ\Xi\rangle$ in (12).

Once $\Xi^{\dagger}$ is determined, the Lagrangian average of a quantity can be computed. For a scalar field $g(x,t)$, in particular, this average is defined by

(Again we keep the notation light by omitting a superscript L for the Lagrangian average.) For an ensemble average this reduces to

We conclude this section with three remarks. First, we give an interpretation for the scalar field $\lambda$. Averaging (35) and comparing with (32) shows that

assuming consistent initial conditions. Letting

(38) is equivalent to

This condition, with $\lambda'$ determined by the condition that $\boldsymbol{u}_{\dagger}$ is divergence free, replaces GLM’s condition $\langle\xi\rangle=0$. Thus, in the polar-factorization framework, the perturbation field has a non-zero mean given by the gradient of $\lambda'$.

The second remark concerns the difference between the bare GLM mean flow map $\bar{\varphi}$ and its volume-preserving approximation $\bar{\varphi}_{\dagger}$. Provided that $\bar{\varphi}$ captures most of the dynamics in the sense that $\Xi$ is close to the identity (or, equivalently, $\xi$ is small), the two maps are close together. We show this explicitly in appendix C where we construct $\bar{\varphi}_{\dagger}$ perturbatively assuming a velocity field of the form

where $\boldsymbol{u}_{0}$ is a fixed background velocity, $\boldsymbol{u}_{1}$ a zero-mean perturbation velocity and $\varepsilon\ll1$. We establish that $\left|\bar{\varphi}(\boldsymbol{a},t)-\bar{\varphi}_{\dagger}(\boldsymbol{a},t)\right|=O(\varepsilon^{2})$ for arbitrarily large $t$. We also consider an alternative volume-preserving mean map proposed by Soward and Roberts [16] and defined perturbatively (see also [17, 25]). This second volume-preserving mean map is also $O(\varepsilon^{2})$-close to $\bar{\varphi}$.

The final remark concerns a third volume-preserving mean map. Instead of approximating the GLM flow map $\bar{\varphi}$ by a volume-preserving flow map as we do for $\bar{\varphi}_{\dagger}$, it is tempting to approximate the GLM mean velocity $\bar{\boldsymbol{u}}$ by a divergence-free velocity field then use the corresponding flow map as volume-preserving mean map. The simplest construction takes the divergence-free velocity field as $P\bar{\boldsymbol{u}}$, the $L_{2}$-projection of $\bar{\boldsymbol{u}}$ on divergence-free vector fields. While this construction is seemingly appealing, we show in appendix C that the resulting flow map does not remain close to $\bar{\varphi}$ for long times, so that it does not properly capture the mean dynamics.

Numerical computation of Lagrangian means

Exponential mean

We now turn to the numerical computation of volume-preserving Lagrangian time averages. We focus on the exponential average

where the parameter $\alpha>0$ is an inverse averaging time scale. The exponential average corresponds to the choice of kernel $k(t)=\alpha e^{-\alpha t}H(t)$, with $H(t)$ the Heaviside function, in the general time average (5). It has the advantage that $\langle f\rangle(t)$ satisfies a differential equation, namely

This makes it possible to compute averages on-the-fly, as part of simulations that deliver $f(t)$ incrementally. Minz et al. [21] exploit this property for the computation of Lagrangian averages and show it generalizes to $k(t)$ given by a sum of exponentials. Baker [26] extends this to time-symmetric kernels by solving differential equations forward and backward in time for the offline computation of Lagrangian averages.

With the average (43), the right-hand side $\langle u \circ\varphi\rangle\circ\bar{\varphi}_{t}^{\dagger}$ of equation (30) governing the evolution of $\lambda$ can be evaluated explicitly. Integration by parts gives

Using (26) and (31) then leads to

Eq. (30) becomes

As in §II B, it is convenient to take the divergence and curl. This gives

These equations should be solved together with the evolution equation (35) for $\Xi_{\dagger}$ and the divergence-free condition (29) for $\bar{u}_{\dagger}$. Eq. (47a) governs the evolution of $\nabla^2\lambda$ and, on inverting the Laplacian, $\lambda$. Eqs. (29) and (35) determine $\bar{u}_{\dagger}$. In an on-the-fly implementation, this system of equations is solved at the same time as dynamical equations for the instantaneous velocity $u$. With $\Xi_{\dagger}$ and $\bar{u}_{\dagger}$ determined, the volume-preserving Lagrangian average of any field can be computed by solving an additional PDE. For a scalar field $g(x,t)$, this takes the form

as follows from differentiating the definition

with respect to $t$.

Volume-preserving Lagrangian averaging based on polar factorisation is not limited to the exponential mean but can be adapted to various types of filters, e.g. those proposed in Refs. [13, 21]. We illustrate this point by developing a volume-preserving version of the Butterworth filter in appendix D.

Application to a two-dimensional incompressible flow

We implement the system above for solutions of the two-dimensional incompressible Euler equations. This system does not admit fast wave motion. We use it solely to compare the impact of standard and volume-preserving Lagrangian averaging on the slow vortical flow. In vorticity–streamfunction formulation, the Euler equations read

where $\zeta$ is the vorticity, $\psi$ the streamfunction and $\nabla^\perp=(-\partial_y,\partial_x)$.

We solve the Euler equations (50) in a $2\pi\times2\pi$ periodic domain using a modified version of the pseudospectral code in Ref. [20]. We use $256^2$ grid points and Fourier modes, and an RK4 integrator with time step $\Delta t=0.005$. For numerical stability, we employ $2/3$ dealiasing and we add hyperviscous dissipation to the vorticity equation by multiplying the Fourier transform $\hat{\zeta}(k,t)$ of the vorticity by $\exp(-\kappa|k|^8\Delta t)$ at each time step. We take $\kappa=2.6\times10^{-14}$.

Together with the Euler equations, we solve the equations required for the computation of Lagrangian averages. For volume preserving averaging, these consist of (35), (47) and (48). To deal with periodic fields, we use $\xi_{\dagger}$ instead of $\Xi_{\dagger}$ and $\lambda'=\lambda-\frac{\alpha}{2}\lVert x\rVert^2$ instead of $\lambda$. We enforce the divergence-free constraint for $\bar{u}_{\dagger}$ by means of a streamfunction, letting

With these changes, (35) becomes

and eq:47 become

At each time step, we update $\xi^\dagger$ and $\lambda'$ using (53) and (54) then compute $\bar{\psi}^\dagger$ by solving (55) using the iteration

where $k = 0, 1, \cdots$ denotes the iterate and $\bar{u}^{\dagger,k} = \nabla^\perp\bar{\psi}^{\dagger,k}$. We focus on the Lagrangian average $\bar{\zeta}^\dagger$ of the vorticity, computed from (49) with $g = \zeta$. We use bilinear interpolation to evaluate $u$ and $g$ at position $\Xi^\dagger(x,t) = x + \xi^\dagger(x,t)$ in (49) and (53). For numerical stability, we find it necessary to add to (53) the same hyperviscous dissipation as in the vorticity equation. For comparison, we also compute the standard GLM average $\bar{\zeta}$ of the vorticity. For this, we solve the bare versions of (49) and (53) with $\bar{u} = \alpha\xi$, as follows from (44) [see also [21]].

The top row of Figure 2 shows the evolution of the vorticity $\zeta$ for an initial condition that consists of two well-separated vortices, specifically

The dynamics are familiar: the vortices orbit the central point $(0,0)$, deform and undergo vortex merger. The second and third rows show the Lagrangian averages $\bar{\zeta}$ and $\bar{\zeta}^\dagger$ of the vorticity obtained with the exponential filter with inverse time scale $\alpha= 0.5$. To mitigate the effect of the phase delay that affects the exponential filter, we offset time by the zero-frequency phase delay $\langle t\rangle= \alpha^{-1} = 2$.

The bare Lagrangian averaging exhibits the expected divergence effect: the vortices in the Lagrangian mean field $\bar{\zeta}$ have much decreased areas compared with those in the instantaneous field $\zeta$. This is successfully corrected by the volume-preserving average. The averaging distorts the vorticity field, primarily because the frequency-dependent phase shift it introduces is only partially compensated by the time offset. Thus, while the slowly evolving filaments at the periphery of the vortex merger are in phase in $\zeta$ and $\bar{\zeta}^\dagger$, the faster evolving cores are not.

The bottom row of Figure 2 gives a finite-amplitude illustration of the difficulties that arise with the simple replacement of the bare mean velocity $\bar{u}$ by its divergence-free projection $P\bar{u}$ discussed for small-amplitude perturbations in appendix C. With this choice, the mean flow map does not stay close to the actual flow map; equivalently, the perturbation map $\Xi$ does not stay close to the identity map. As a result, the mean vorticity field develops small-scale features unrelated to the instantaneous dynamics.

A complementary view of the divergence effect and its correction by volume-preserving averaging is offered in Figure 3. The figure shows the probability density function (pdf) of $\zeta$, $\bar{\zeta}$ and $\bar{\zeta}^\dagger$ at $t = 20$. The pdf of $\zeta$ is time independent up to numerical effects (truncation error and dissipation). The pdfs of $\bar{\zeta}$ and $\bar{\zeta}^\dagger$ are also time independent after an initial adjustment phase. The depletion of high values of $\bar{\zeta}$ reflects the decrease in area of the vortices noted in Figure 2. The pdf of $\bar{\zeta}^\dagger$, in contrast, is almost identical to that of $\zeta$, confirming the effectiveness of the volume-preservation algorithm.

Application to a rotating shallow-water flow

In this section, we apply Lagrangian averaging to a two-time-scale flow, namely a rapidly rotating shallow-water flow that combines slow turbulent vortical motion with a fast mode-1 Poincaré wave. Such a flow is used in Refs. [13, 19, 21] to demonstrate the capabilities of different PDE-based averaging methods.

At first sight, the divergence effect that we address in the present paper does concern the shallow-water model. In this model, the (two-dimensional, depth-averaged) velocity field is divergent because of the free surface, so there appears to be no reason to insist that the mean flow should preserve areas. However, in the rapidly rotating, small-Rossby-number regime that we consider, the vortical part of the flow is approximately in geostrophic balance. This implies that the corresponding velocity is approximately divergence free. More precisely, the Lagrangian mean velocity is approximately geostrophically balanced and divergence free [see [23], and references therein]. For sufficiently small Rossby number, we expect the divergence effect associated with the definition of the Lagrangian mean flow to dominate the small physical divergence associated with finite Rossby number. It is therefore sensible to impose exact area preservation in order to eliminate the divergence effect.

We carry out a simulation using the rotating shallow water model, as given in [21]. The non-dimensional momentum and continuity equations are [27]

where $Ro = U/(fL)$ is the Rossby number and $Fr = U/\sqrt{gH}$ is the Froude number. Here $U$, $L$ and $H$ are the reference velocity, length and depth used for the non-dimensionalization.

We initialise the model with the same initial condition as in [20] and [21], superimposing a geostrophic turbulent flow (obtained by prior solution of an incompressible two-dimensional Navier–Stokes model for the velocity and geostrophic balance for the height field) and a right-travelling mode-1 Poincaré wave. We take the root-mean square velocity of the geostrophic flow as reference velocity $U$ and the length scale of the first Fourier mode as reference length $L$ so that the doubly periodic domain is $[0, 2\pi]^2$. The Rossby and Froude number of the geostrophic flow are $Ro = 0.1$ and $Fr = 0.5$. With these parameters, the frequency $\omega= (Ro^{-2} + Fr^{-2})^{1/2}$ of the mode-1 Poincaré wave is $\omega= 10.2$ corresponding to a period of 0.62 time units. The wave fields in the absence of flow are given by

At the initial time $t = 0$, we add these to the geostrophic flow field. We take the amplitude $a = -1$ so that the maximum wave velocity is as large as the geostrophic flow root-mean-square velocity.

To obtain the volume preserving Lagrangian mean of a scalar field, we solve (48) along with the dynamical equations (56), equation (52) for $\xi^{\dagger}$ and equations (53) for $\lambda'$ and $\psi^{\dagger}$. We use a pseudospectral implementation, identical to that of §III B for (52) and (53). For numerical stability, we apply the same dissipation scheme to (56a) and (52) as in §III B.

Figure 4 shows the evolution of the vorticity $\zeta$ and of its bare and volume-preserving averages $\bar{\zeta}$ and $\bar{\zeta}^{\dagger}$. The instantaneous vorticity is dominated by the large-scale wave signal, but vortical structures are visible as smaller scale modulations. Much of the wave signal is eliminated by the averaging, leaving similar-looking fields $\bar{\zeta}$ and $\bar{\zeta}^{\dagger}$. The vortices in $\bar{\zeta}$ are smaller than those in $\bar{\zeta}^{\dagger}$ and, though this is less obvious to discern, in $\zeta$. Vorticity (as opposed to potential vorticity) is not materially conserved in shallow water, so the area of vortices change in time. Nonetheless we argue that the area changes in $\bar{\zeta}$ are the result of the divergence effect rather than of a physical divergence.

To support this, we show in figure 5 the divergences $\nabla\cdot\langle\mathbf{u}\rangle$ and $\nabla\cdot\bar{\mathbf{u}}$ of the Eulerian and bare Lagrangian mean velocities. The bare Lagrangian mean $\bar{\mathbf{u}}$ for the exponential mean is computed from the relation $\bar{\mathbf{u}} = \alpha(\mathbf{\Xi} - \mathbf{x})$ (which is shown in [21] or can be deduced from (46) by setting $\lambda= |\mathbf{x}|^2/2$). Both $\nabla\cdot\langle\mathbf{u}\rangle$ and $\nabla\cdot\bar{\mathbf{u}}$ are dominated by a wave signal that results from the limitation of exponential averaging as a low-pass filter. The slow signal, correlated with the vortical structure in figure 4, is much larger in $\nabla\cdot\bar{\mathbf{u}}$, suggesting it results from the averaging rather than from a physical process. We also show in figure 5 (right panel) the Jacobian field $|\partial\mathbf{\Xi}(\mathbf{x},t)/\partial\mathbf{x}|$. Since

this Jacobian is the ratio of the height changes associated with the compressibility of $\bar{\varphi}$ and $\varphi$, respectively. The compressibility of $\varphi$ is weak because $\nabla\cdot u$ is dominated by an oscillatory wave component that has little time-integrated effect. The large values in $\left|\partial\Xi(x,t)/\partial x\right|$, up to 2, are dominated by the compressibility of $\bar{\varphi}$ induced by the divergence effect.

Figure 6 shows the first components $u$, $\bar{u}$ and $\bar{u}^{\dagger}$ of the velocity and of its Lagrangian means at fixed $y=0.24$ as functions of $t$ and $x$. The instantaneous $u$ (top row) is dominated by the fast Poincaré wave. The bare Lagrangian mean $\bar{u}$ primarily captures the slow dynamics of the simulation, but a substantial wave signal remains, reflecting the limitation of the exponential mean as a low-pass filter. Remarkably, the wave signal is much reduced in the volume-preserving version $\bar{u}^{\dagger}$ (third row) obtained from eq:47. We attribute this improvement to the fact that the Poincaré wave is divergent and partially filtered from $\bar{u}^{\dagger}$ by the divergence-free constraint $\nabla\cdot u^{\dagger}=0$. We find the same improvement between the bare and volume-preserving Lagrangian averages when the exponential filter is replaced by the 2nd-order Butterworth filter (fourth and fifth rows of Figure 6). The benefit of volume preservation combines with the better low-pass filtering of the Butterworth filter to produce a Lagrangian mean velocity $\bar{u}^{\dagger}$ that is almost perfectly free of fast oscillations.

Discussion

This paper develops and implements a method for the computation of Lagrangian averages referred to mean positions that are related to the initial positions of particles by a volume-preserving mean flow map $\bar{\varphi}^{\dagger}$. Equivalently, the Lagrangian mean velocity $\bar{u}^{\dagger} = \partial_t\bar{\varphi}^{\dagger} \circ\bar{\varphi}^{\dagger-1}$ is divergence free: $\nabla\cdot\bar{u}^{\dagger} = 0$. For incompressible fluids, this also means that the perturbation map $\Xi^{\dagger}$ between mean and actual positions preserves volume.

The volume-preservation properties are in contrast with the standard GLM for which the mean and perturbation maps $\bar{\varphi}$ and $\Xi$ do not preserve volume and the Lagrangian mean velocity is non-divergent, $\nabla\cdot\bar{u} \ne0$. Our construction is nonetheless anchored to standard GLM since $\bar{\varphi}^{\dagger}$ is chosen as the volume-preserving map closest to $\bar{\varphi}$ or, equivalently, as the volume-preserving factor in the polar factorization of $\bar{\varphi}$. Note however that $\bar{\varphi}^{\dagger}$ can be alternatively defined as the Fréchet mean of an ensemble of flow maps, without explicit reference to the GLM mean map (see (24)).

The polar-factorization characterization we use relies on the (arbitrary) choice of the $L_2$-norm to measure the closeness of $\bar{\varphi}^{\dagger}$ to $\bar{\varphi}$. In principle, other $L_p$-norms could be chosen, but the corresponding optimal-transport theory is considerably more complicated. Alternatively, as proposed in [17], it is possible to replace the pointwise $L_2$-norm by a more intrinsic measure of distance between maps, namely the distance that arises naturally in the interpretation of flow maps as points in a Lie group with right-invariant metric [e.g. [28]]. Turning the resulting abstract definition of the mean map (termed “geodesic mean” in [17]) into a practical algorithm is a formidable challenge, however. The volume-preserving mean of Soward and Roberts [16] (termed “glm”, see [17]), which is defined perturbatively assuming small-amplitude perturbations, is similarly challenging. With this in mind, we contend that the polar-factorization definition of a volume-preserving mean map offers the most convenient way of overcoming the divergence effect of GLM in numerical Lagrangian averaging.

Our simulation results suggest that the divergence effects can strongly distort the slow dynamics. This may not have been fully appreciated in the past because most studies of GLM assume small-amplitude perturbations and track the slow dynamics over relatively short times. Over multiple eddy times scales, however, the divergence effect leads to a significant contraction of vortices and, conversely, significant expansion of hyperbolic regions. This occurs even though the divergence of the Lagrangian mean velocity is weak because the process of volume or area change is multiplicative. The distortion of the slow dynamics is problematic when Lagrangian averaging is used to separate fast waves from slow flows by substracting full and averaged fields [e.g. [13]]. It would be worthwhile investigating how the adoption of the volume-preserving mean flow of this paper improves this separation.

This paper focuses entirely on the kinematics of GLM and its volume-preserving variant. Key advantages of GLM over Eulerian averaging lie however in the structure of its averaged dynamical equations. An important question is therefore whether this structure persists for volume-preserving Lagrangian averaging. The answer is yes. As Gilbert and Vanneste [17, 25] emphasise, much of the structure of the GLM equations is independent of the specific choice made for the Lagrangian mean map. The Lagrangian mean momentum equation, in particular, is only affected through changes in the specific form of the pseudomomentum, with extra terms arising because $\langle\zeta\rangle\ne0$. Whether the associated complications outweigh the benefits brought about by volume preservation may be a matter ot taste.

Supplementary movies are available at https://abhijeet-minz.github.io/animations.html

The code developed for the numerical applications in §III B and §III C and is openly available [29].

We are grateful to H. A. Kafiabad and C. Maitland-Davies for useful discussions.

The authors report no conflict of interest.

Rotating particle

Consider a two-dimensional fluid in rigid rotation with angular velocity $\Omega(t)$. The position of a particle, $z(t) = x(t) + \mathrm{i}y(t)$ in complex notation, satisfies

and is given by

where

The initial position $c=a+\mathrm{i}b$ can be interpreted as a label for the particle. The corresponding Lagrangian mean position is

We focus on time averages of the form (5). With the additional assumption $k(t)\ge0$, the triangle inequality implies that $\lvert\chi\rvert\le1$, hence

This shows that averaging contracts trajectories. The Lagrangian mean velocity

is divergent,

unless $\dot{\theta}=\mathrm{const.}$

In this paper we propose a definition of the mean flow map that replaces the straightforward coordinate-wise mean of GLM by the closest volume-preserving flow map. For rigid-body rotation, this amounts to replacing $\bar{z}(t)$ in (65) by the nearest point on the circle of radius $\lvert z(t)\rvert=\lvert c\rvert$, that is, by

This represents a rigid-body rotation. The velocity

is divergence-free.

Polar factorization

We give a formal argument for the polar factorization (26). First observe that

where the constant is independent of $\psi$ because

using that $\psi$ preserves volume. We next consider the maps $\psi$ as near-identity, volume-preserving perturbations of the extremum $\bar{\phi}^{\dagger}$. Near-identity perturbations can be expressed in terms of a generating vector field $w$ as

The divergence-free condition $\nabla\cdot w=0$ ensures that these perturbations preserve volume. With $\psi=e^{w\cdot\nabla}\circ\bar{\phi}^{\dagger}$, this gives

With this expression and (71), the condition that $\psi=\bar{\phi}^{\dagger}$ minimises the left-hand side of (71) reads

on changing integration variable from $a$ to $x = \bar{\varphi}_{\dagger}^{-1}(a,t)$ and omitting $O(\lVert w\rVert^3)$ terms. For the inequality to hold, the first integral on the last line of (75) must vanish for all divergence-free $w$. This implies that

for some $\lambda$, as follows from the orthogonality of the Helmholtz decomposition [e.g. [30]]. Rearranging gives the polar factorization (26). Introducing (76) into the second integral and integrating by parts gives

i.e. convexity of $\lambda$.

Small-amplitude perturbations

We now obtain an approximation to the volume-preserving Lagrangian mean velocity $\bar{u}_{\dagger}$ for a velocity field of the form (41). The velocity is divergence free such that $\nabla\cdot u_0 = \nabla\cdot u_1 = 0$. We use an ensemble average, with $\langle u_1\rangle= 0$. We solve (29), (32) and (35) perturbatively. We expand

anticipating some $O(1)$ and $O(\varepsilon)$ terms. Introducing this into (35) gives

at $O(\varepsilon)$. Here we introduce the notation

for the Lie derivative (or commutator) of vector fields. We note that the solution $\xi_1$ satisfies $\langle\xi_1\rangle= 0$ and $\nabla\cdot\xi_1 = 0$. For divergence-free vector fields like $u_0$ and $\xi_1$, the Lie derivative takes the alternative form

making it plain it is also divergence free. Averaging the $O(\varepsilon^2)$ terms in (35) yields

where we use that $\langle\xi\rangle_1 = 0$. The middle term on the right-hand side of (84) has the component expression

with summation over repeated indices implied. This term can be rewritten by noting that

This identity can established in components starting with the first term on the right-hand side:

using (81) in the form $(\partial_t + u_0^j\partial_j)\xi_1^k = u_1^k + \xi_1^j\partial_j u_0^k$. Rearranging yields (86).

Introducing (86) into (84) gives

Now, according to (40), $\langle\xi_2\rangle$ is a gradient, $\langle\xi_2\rangle=\nabla\lambda_2$, with $\lambda_2$ determined by the condition $\nabla\cdot\bar{u}_{\dagger2}=0$. Since $\nabla\cdot\langle\mathcal{L}_{\xi_1}u_1\rangle=0$, this implies

where $P$ denotes the (Leray) projection on divergence-free vector fields. This reduces (88) to

and recovers the result of the ‘optimal transport’ formulation of [17].

We can now compare the volume-preserving Lagrangian mean velocity $\bar{u}^{\dagger}=u_0+\varepsilon^2\bar{u}_{\dagger2}$ resulting form the polar factorization (26) with other Lagrangian mean velocities. The standard GLM Lagrangian mean velocity is $\bar{u}=u_0+\varepsilon^2\bar{u}_2$ with

It is obtained by setting $\langle\xi_2\rangle=0$ in (88). Comparing with (90) shows that the polar factorization construction introduces the projection $P$ in the last term. The glm construction of Soward and Roberts [16], reducing to the solenoidal Stokes drift of Vanneste and Young [18] for $u_0=0$, is

An important point is that the three Lagrangian mean velocities (90), (91) and (92) differ by terms that are derivatives along the leading order flow, $\partial_t+\mathcal{L}_{u_0}$. This implies that the mean trajectories corresponding to the three mean velocities are $O(\varepsilon^2)$ close for long times. More precisely, we have that

where $\bar{\phi}^{\mathrm{glm}}$ is the flow map associated with $u_0+\varepsilon^2\bar{u}^{\mathrm{glm}}_2$ and $f=\frac{1}{2}P\langle\xi_1\cdot\nabla\xi_1\rangle$.

We verify (93) by taking the time derivative of the right-hand side to obtain

This matches the time derivative $\partial_t\bar{\phi}^{\dagger}(a,t)=\bar{u}^{\dagger}(\bar{\phi}^{\dagger}(a,t),t)$ of the left-hand side.

A similar computation shows that

with $f=\frac{1}{2}\langle\xi_1\cdot\nabla\xi_1\rangle$.

We conclude from the closeness of $\bar{\phi}$, $\bar{\phi}^{\dagger}$ and $\bar{\phi}^{\mathrm{glm}}$ shown by (93) and (95) that they represent the mean motion equally well.

We contrast this closeness with the behaviour of the trajectories of a superficially appealing divergence-free alternative to $\bar{u}$, namely its projection $P\bar{u}=u_0+\varepsilon^2P\bar{u}_2$ on divergence-free vector fields. Since

and $P$ and $\mathcal{L}_{u_0}$ do not commute, the difference between $P\bar{u}$ and $\bar{u}$ is not a derivative along the leading order flow. As a result, the trajectories of $P\bar{u}$ do not remain close to $\bar{\phi}$; they do not capture the mean motion and lead to the uncontrolled growth of the displacement field $\xi$.

Volume-preserving Lagrangian averaging with a Butterworth filter

Volume-preserving Lagrangian averaging is not limited to the exponential filter of §III. It can be extended to a broad class of filters including the top-hat filter of Ref. [20] or the 2nd-order Butterworth filter of Ref. [21]. The key step is to identify the bare mean flow map $\langle\varphi\rangle=\bar{\varphi}$ and substitute the polar decomposition $\nabla\lambda\circ\bar{\varphi}_{\dagger}^{-1}$.

To illustrate this, we apply the construction to the 2nd-order Butterworth filter, an efficient low-pass filter in the class of sum-of-exponential filters considered by Minz et al. [21]. For this filter, the mean flow map $\bar{\varphi}$ can be obtained by solving the system of equations

where $\widetilde{\varphi}$ is an auxiliary field and the matrix

is positive definite. Defining $\widetilde{\Xi}=x+\widetilde{\xi}$ by $\widetilde{\varphi}=\widetilde{\Xi}\circ\bar{\varphi}$, this can be rewritten as

where $\widetilde{u}=\alpha\bigl(\xi-(\sqrt{2}-1)\widetilde{\xi}\bigr)$ and $\bar{u}=\alpha\widetilde{\xi}$.

The auxiliary velocity $\widetilde{u}$ is defined by $\widetilde{u}=\partial_t\widetilde{\varphi}\circ\bar{\varphi}^{-1}$ (see [21]).

We now obtain a volume-preserving version of eq:D3. We substitute $\bar{\varphi}=\nabla\lambda\circ\bar{\varphi}_{\dagger}$ into (97) and compose with $\bar{\varphi}_{\dagger}^{-1}$ to find

where $\widetilde{u}_{\dagger}=\partial_t\widetilde{\varphi}\circ\bar{\varphi}_{\dagger}^{-1}$, $\widetilde{\Xi}_{\dagger}=\widetilde{\varphi}\circ\bar{\varphi}_{\dagger}^{-1}$ and $\Xi_{\dagger}=\varphi\circ\bar{\varphi}_{\dagger}^{-1}$. With $\widetilde{\Xi}_{\dagger}=x+\widetilde{\xi}_{\dagger}$ and $\lambda=\lvert x\rvert^2/2+\lambda'$, the first row in (101) becomes

where $\widetilde{\Xi}_{\dagger}=x+\widetilde{\xi}_{\dagger}$ and $\Xi_{\dagger}=x+\xi_{\dagger}$. The second row of (101) is

Together with the divergence-free condition $\nabla\cdot\bar{u}_{\dagger}=0$, (102) and (103) determine $\widetilde{u}_{\dagger}$, $\bar{u}_{\dagger}$ and $\lambda'$ assuming $\widetilde{\xi}_{\dagger}$ and $\xi_{\dagger}$ are known. Evolution equations for $\widetilde{\xi}_{\dagger}$ and $\xi_{\dagger}$ are obtained by time differentiation of $\widetilde{\varphi}=\widetilde{\Xi}_{\dagger}\circ\bar{\varphi}_{\dagger}$ and $\varphi=\Xi_{\dagger}\circ\bar{\varphi}_{\dagger}$ and are given by (99) with all fields decorated by $\dagger$.

In two dimensions, we can take $\bar{u}_{\dagger}=\nabla^{\perp}\bar{\psi}_{\dagger}$. The divergence and curl of (103) then give the system

analogous to (53).

References