*Review* **Potential Fields in Fluid Mechanics: A Review of Two Classical Approaches and Related Recent Advances**

#### **Markus Scholle 1,\*,†, Florian Marner 2,† and Philip H. Gaskell 3,†**


Received: 24 March 2020; Accepted: 24 April 2020; Published: 27 April 2020

**Abstract:** The use of potential fields in fluid dynamics is retraced, ranging from classical potential theory to recent developments in this evergreen research field. The focus is centred on two major approaches and their advancements: (i) the Clebsch transformation and (ii) the classical complex variable method utilising Airy's stress function, which can be generalised to a first integral methodology based on the introduction of a tensor potential and parallels drawn with Maxwell's theory. Basic questions relating to the existence and gauge freedoms of the potential fields and the satisfaction of the boundary conditions required for closure are addressed; with respect to (i), the properties of self-adjointness and Galilean invariance are of particular interest. The application and use of both approaches is explored through the solution of four purposely selected problems; three of which are tractable analytically, the fourth requiring a numerical solution. In all cases, the results obtained are found to be in excellent agreement with corresponding solutions available in the open literature.

**Keywords:** potential fields; Clebsch variables; Airy's stress function; Goursat functions; Galilean invariance; variational principles; boundary conditions; film flows

**PACS:** 47.10.-g

**MSC:** 76M99

#### **1. Introduction**

In various branches of physics, potentials continue to be used as additional auxiliary fields for the advantageous reformulation of one or more governing equations. Classical electrodynamics serves as a prime illustrative example of their use [1], enabling Maxwell equations in their commonly expressed form, comprised of two scalar and two vector equations for the two observables—electric field *E* and magnetic flux density *B*—to be formulated differently in terms of the derivatives of a scalar potential *ϕ* and a vector potential *A* , such that *<sup>E</sup>* <sup>=</sup> −∇*<sup>ϕ</sup>* <sup>−</sup> *<sup>∂</sup>tA* and *<sup>B</sup>* <sup>=</sup> ∇ × *<sup>A</sup>* , respectively; in which case, two of the four Maxwell equations are fulfilled identically while the other two form a self-adjoint pair—i.e., can be obtained from a variational principle. Moreover, followed by proper gauging of the potential fields, a fully decoupled form of Maxwell's equations in terms of two inhomogeneous d'Alembert equations is obtained. From a mathematical viewpoint, amendments to the original equations comprise:


Despite the above benefits, a historical debate reaching back to Heaviside [2,3] has surrounded the Maxwell potentials and their rationale, which were not free of polemics [4–6]. A major objection against the use of potential fields has surrounded their mathematical non-uniqueness, implying that they are not observable by physical measurements. In the case of the Maxwell potentials, the transformation

$$
\vec{A} \to \vec{A}' = \vec{A} + \nabla \chi \tag{1}
$$

$$
\varphi \to \varphi' = \varphi - \partial\_l \chi \tag{2}
$$

in terms of an arbitrary scalar field *χ* defines an alternative set of potentials related identically to the observable fields *E* and *B*. The non-uniqueness of the original Maxwell potentials does not reflect a weak point in the theory; on the contrary, it defines a symmetry that has proven formative in the usage of potential fields in modern theoretical physics as well as ground-breaking with respect to the subsequent development of gauge theory [7]. Each of which confirm Maxwell's theory to be a significant reference point for other field theories, an obvious example being the potential representation of the Weyl tensor in general relativity [8,9].

Subsequently, based on the publications of Ehrenberg and Siday [10], Aharonov and Bohm [11], the so-called *Aharonov–Bohm solenoid effect*—which takes place when the wave function of a charged particle passes around a long solenoid and experiences a phase shift as a result of the enclosed magnetic field, despite the magnetic field being negligible outside the solenoid—gave rise to a second debate concerning the physical meaning of the vector Potential *A* . This effect is frequently misunderstood: it does not allow for a point-wise "measurement" of *A* since only the integral magnetic flux can be determined from the phase shift, while the gauge transformation in Equation (1) is still a symmetry of the Schrödinger–Maxwell theory predicting the effect properly; thus any gradient field ∇*χ* can be added to the vector potential *A* leading to the same experimental result. Nevertheless, the question still remains an open one as to what physical role the vector potential plays and controversy surrounding it persists; see, for example, [12–18].

In summary, the use of potentials is motivated from both a mathematical and physical viewpoint: mathematically the original equations can be manipulated in beneficial ways, while physically new insight concerning the structure of the theory is provided via the reformulated equation set. Both of these aspects are considered subsequently.

The particular aim henceforth is to demonstrate the application of potential methods to fluid mechanics. As is well known in the case of an incompressible inviscid irrotational fluid flow, by expressing the velocity field *v* in terms of a scalar potential *ϕ*:

$$
\vec{v} = \nabla \varphi \tag{3}
$$

the associated equations of motion are reduced to Bernoulli's equation, resulting as a *first integral* of Euler's equations together with a Laplace equation

$$
\nabla^2 \varphi = 0 \tag{4}
$$

for the potential via the continuity equation. Since Bernoulli's equation is essentially the conditional equation for the pressure field, only Equation (4) remains to be solved, elevating potential flow theory to the status of an essential topic in standard fluid dynamics text books [19–22]. Despite the obvious advantage of making various flow problems more tractable, the approach is restricted to inviscid and irrotational barotropic flows. In order to extend the methodology to encompass a wider range of flow problems, generalisations of and alternatives to Equation (3) have emerged. Of these, the following two major strands only are explored in the present work:


Historically, both approaches have been subject to limitations on their usage: the Clebsch transformation originated for the case of inviscid flow (Re → ∞) and the complex variable method for that of 2D Stokes flow (Re → 0). Recent advancements have lifted these restrictions. In this review, the origin and evolution of both methods is retraced and potential future developments highlighted. Section 2 considers the Clebsch transformation, starting from its origins in inviscid barotropic flow theory, Section 2.1. After commenting on the global existence of the Clebsch variables, Section 2.2, it is demonstrated via a rigorous symmetry analysis, Section 2.2, that the Clebsch representation of the velocity is a natural outcome of Galilean invariance via Noether's theorem. An extended Clebsch transformation for viscous flow is then presented in Section 2.3 and followed, Section 2.4, by applying it to the problem of viscous stagnation flow. The complex variable method and its progression to a tensor potential approach is outlined in Section 3, beginning in Section 3.1 with the use of Airy's stress function with respect to steady-state equilibrium conditions for an arbitrary continuum in general and for Stokes's flow in particular. In Section 3.2, the approach is generalised to encompass the full unsteady 2D-NS equations, utilising a complex potential. A serendipitous benefit is that of enabling the integration of the dynamic boundary condition along a free surface, or interface, and its reduction to a standard form, as shown in Section 3.3, as subsequently utilised in Section 3.4 for the numerical solution of a free surface film flow problem. Latterly, a particularly noteworthy advance, Section 3.5, has been that of overcoming the 2D restriction. This has been achieved by rearranging the complex equations in a tensor form leading to a potential-based first integral of the full 3D-NS equations with a seamless extension to an analogous form of the 4D relativistic energy momentum equations. Concluding remarks are provided in Section 4.

Finally, it would be remiss and incomplete not to point out that in the field of fluid mechanics other approaches to the use of potentials for solving the equations of motion exist that have not been considered in the present text. In this sense and as instructive examples, the reader is referred to the work of Papkovich and Neuber [40], Lee et al. [41], Greengard and Jiang [42].

#### **2. Clebsch Transformation Approach**

#### *2.1. The Clebsch Transformation for Inviscid Flows*

For inviscid flow, Clebsch [19,20,23] proposed a non-standard potential representation for the velocity field, in terms of the so-called Clebsch variables *ϕ*, *α*, *β*, of the form:

$$
\vec{u} = \nabla \varphi + \mathfrak{a} \nabla \beta \,. \tag{5}
$$

From a mathematical viewpoint, the potential representation in Equation (5) is a decomposition of the velocity field into a curl-free part ∇*ϕ* and a helicity-free part *α*∇*β*. Schoenberg [29] has shown that the above decomposition is not unique; by applying the gauge transformation:

$$\begin{aligned} \varphi &\longrightarrow & \varphi' &= \varphi + f(\alpha, \beta, t) \\ \alpha &\longrightarrow & \alpha' &= g(\alpha, \beta, t) \\ \beta &\longrightarrow & \beta' &= h(\alpha, \beta, t) \end{aligned} \tag{6}$$

an equivalent set of Clebsch variables *ϕ* , *α* , *β* is given if and only if the functions *f* , *g*, *h* fulfil the following two PDEs:

$$\frac{\partial f}{\partial \beta} + g \frac{\partial h}{\partial \beta} = \quad \text{a.} \tag{7}$$

$$\begin{aligned} \frac{\partial f}{\partial \alpha} + \mathbf{g} \frac{\partial h}{\partial \alpha} &= \quad 0 \end{aligned} \tag{8}$$

By applying Equation (5) to Euler's equations for inviscid flows, one obtains:

$$\vec{0} = \frac{\text{D}\vec{u}}{\text{D}t} + \nabla[P + \mathcal{U}] = \nabla\left[\frac{\partial\rho}{\partial t} + a\frac{\partial\beta}{\partial t} + \frac{\vec{u}^2}{2} + P\left(\boldsymbol{\varrho}\right) + \mathcal{U}\right] + \frac{\text{Da}}{\text{D}t}\nabla\beta - \frac{\text{D}\beta}{\text{D}t}\nabla a \,, \tag{9}$$

with the pressure function *P* () = ´ −1d*p* and *U* the specific potential energy of the external force. The operator D/D*t* = *∂*/*∂t* +*u* · ∇ is the material time derivative. Being basically of the form:

$$
\nabla \left[ \cdot \cdot \cdot \right] + \left[ \cdot \cdot \cdot \right] \nabla \mathfrak{a} + \left[ \cdot \cdot \cdot \right] \nabla \beta = \mathfrak{d} \,, \tag{10}
$$

this vector equation can be decomposed according to:

$$\frac{\partial \mathcal{Q}}{\partial t} + a \frac{\partial \mathcal{B}}{\partial t} + \frac{i\vec{l}^2}{2} + P + \mathcal{U} \quad = \quad F(\mathfrak{a}, \mathfrak{F}, \mathfrak{t}), \tag{11}$$

$$\frac{\text{D}\mathfrak{a}}{\text{D}t} \quad = \quad -\frac{\partial F}{\partial \beta} \, \tag{12}$$

$$\frac{\text{D}\emptyset}{\text{D}t}^{\text{D}} = \begin{array}{c} \frac{\text{\partial F}}{\partial \alpha'} \\ \end{array} \tag{13}$$

containing an unknown function *F*(*α*, *β*, *t*); which, making use of the gauge transformation of Equation (6), leads to *F* → 0. The above three scalar field equations represent a first integral of Euler's equations and are self-adjoint. However, their most intriguing feature is that the vorticity:

$$
\vec{\omega} = \frac{1}{2}\nabla \times \vec{u} = \frac{1}{2}\nabla \vec{u} \times \nabla \beta,\tag{14}
$$

is given by the two scalar fields *α*, *β*, only. Hence, the vortex dynamics is reduced to and conveniently captured by the two transport Equations (12) and (13). This beneficial feature has been exploited by, for example, Prakash et al. [43], who utilised the Clebsch transformation as a generalisation of classical potential theory for the simulation of bubble dynamics.

Another beneficial feature of the transformed field equations is their *self-adjointness*: they result from a variational principle:

$$\delta \int\_{t\_1}^{t\_2} \iiint\_V \ell \left( \wp, a, \eth \wp, \eth \wp, \ddagger, t \right) \mathrm{d}V \mathrm{d}t = 0,\tag{15}$$

with respect to free and independent variation of the four fields , *ϕ*, *α*, *β* and their first order spatial and temporal derivatives for fixed values at initial and final time, *t*1,2, with the Lagrangian given as [25]:

$$\ell\left(\boldsymbol{\varrho},\boldsymbol{u},\partial\boldsymbol{\varrho},\partial\boldsymbol{\beta},\vec{\mathbf{x}},t\right) = -\varrho\left[\frac{\partial\boldsymbol{\varrho}}{\partial t} + \boldsymbol{u}\frac{\partial\boldsymbol{\beta}}{\partial t} + \frac{1}{2}\left(\nabla\boldsymbol{\varrho} + \boldsymbol{u}\nabla\boldsymbol{\beta}\right)^2 + \boldsymbol{e}\left(\boldsymbol{\varrho}\right) + \boldsymbol{U}\left(\vec{\mathbf{x}},t\right)\right],\tag{16}$$

for the specific elastic energy *e*() fulfilling *e*() + *e* () = *P* (). For generalisations of the above variational principle toward thermal degrees of freedom, reference is made to [24,44–46].

#### *2.2. A Note on the Global Existence of the Clebsch Variables*

When introducing potentials as auxiliary fields, their existence has to be clarified first. Referring to the classical example *u* = ∇*ϕ* from potential theory corresponding to the special case of the Clebsch representation in Equation (5) with *α* = *β* = 0, it is obvious that it applies only to vortex-free fields, i.e., ∇ × *<sup>u</sup>* <sup>=</sup> 0. However, according to Equation (14), velocity fields with non-vanishing vorticity (Equation ((14)) can be considered, but not arbitrary ones, as demonstrated in the following, Moreau [47], Moffatt [48] having identified the *helicity*

$$H = \iiint\limits\_{V} 2\vec{u} \cdot \vec{\omega} \,\mathrm{d}V \tag{17}$$

as a decisive quantity in the case of inviscid fluid flow. Under the assumption that the potential *ϕ* is continuously differentiable and single-valued then by expressing, via Equations (5) and (14), the velocity and the vorticity in terms of the Clebsch variables and their gradients, the helicity density can be re-written as:

$$\begin{aligned} 2\vec{u} \cdot \vec{\omega} &= \left[\nabla\varrho + \mathfrak{a}\nabla\beta\right] \cdot \left(\nabla\mathfrak{a} \times \nabla\beta\right) = \nabla\varrho \cdot \left(\nabla\mathfrak{a} \times \nabla\beta\right) \cdot \vec{\omega} \\ &= \nabla\varrho \cdot \vec{\omega} = 2\nabla \cdot \left(\varrho\vec{\omega}\right) - 2\varrho\underbrace{\nabla\cdot\vec{\omega}}\_{\vec{0}} \end{aligned}$$

implying, utilising Gauss's theorem with *∂V* denoting the surface of the flow domain and*n* the normal vector, the global helicity to be:

$$H = 2 \iiint\limits\_{V} \nabla \cdot (q \vec{\omega}) \,\mathrm{d}V = \iiint\limits\_{\partial V} q \vec{\omega} \cdot \vec{n} \mathrm{d}A \,. \tag{18}$$

The issue raised by the above formula is the occurrence of the potential *ϕ* as a non-observable in the sense that it can, according to Equation (6), be replaced by a re-gauged potential *ϕ* = *ϕ* + *f* (*α*, *β*), leading to a corresponding helicity given by:

$$H' = H + 2\iiint\limits\_{\partial V} f\left(\mathfrak{a}, \mathfrak{k}\right) \vec{\omega} \cdot \vec{n} \mathrm{d}A\,. \tag{19}$$

Since, due to its definition (Equation ((17)), the helicity is an observable and therefore a gauge-invariant quantity, i.e., *H* = *H*, the integral on the right hand of Equation (19) has to vanish for any choice of the function *f* (*α*, *β*), implying that *ω* · *n* = 0 along the entire boundary *∂V* of the flow domain. Finally, the latter implies also that the helicity vanishes. As a consequence, the classical Clebsch transformation with continuously differentiable and single-valued potentials only applies to flows the total helicity of which is zero.

There are two different approaches to overcoming this restriction: the first utilises a multi-valued potential *ϕ*, as demonstrated by Yahalom [49], Yahalom and Lynden-Bell [50]; the second is based on the use of multiple pairs of variables, such as *u* = ∇*ϕ* + *α*1∇*β*<sup>1</sup> + *α*2∇*β*<sup>2</sup> + ··· , in the sense that the number of pairs has to be chosen adequately depending on the topological features of the respective individual flow problem—such flows with closed vortex lines that form linked rings or ones with isolated points of zero vorticity are discussed, for example, by Balkovsky [51], Yoshida [52]. More recent work on this topic is presented by Ohkitani and Constantin [53], Cartes et al. [54], Ohkitani [55]. However, independently of the question of how many pairs of Clebsch variables are useful for representing the topology of a flow, there is a minimum number of two pairs for which global existence can be granted, as demonstrated in Appendix A.

Subsequently, for convenience, attention is paid to the classical form in Equation (5) only on the understanding that an extension via more pairs of variables is possible if required.

#### Derivation of a Clebsch-Like Form by Galilean Invariance and Self-Adjointness

The question of global existence discussed briefly above, Section 2.2, has become a research topic to which many research groups have contributed over several decades, while the use of the Clebsch variables remained restricted to the realm of inviscid flows. In contrast, the focus here is to provide a rationale for the use of Clebsch variables for arbitrary continuous physical systems, the dynamics of which are assumed to be given by a variational principle, following an in-depth analysis of the underlying Galilean group and its consequences for the resulting balances of mass and momentum.

If a system is *physically closed*, i.e., isolated from the surrounding, its equations of motion are invariant with respect to the following four universal symmetry transformations of the Galilean group, corresponding to homogeneity of time and space, isotropy of space and equality of all inertial frames:

#### **time translations:**

**space translations:**

*t* → *t* = *t* + *τ* (20)


**Galilei boosts:**

$$
\vec{\mathbf{x}} \to \vec{\mathbf{x}}' = \vec{\mathbf{x}} - \vec{u}\_0 t \tag{23}
$$

Here, the scalar *τ*, the two vectors*s* and *u*<sup>0</sup> and the unitary matrix R fulfilling R*T*R = 1 and det R = 1 are constants. Via Formulae (20)–(23) the four symmetries are obviously well-defined for discrete systems; for instance, systems of point masses in Newtonian mechanics. For continuous systems, the situation is essentially different: in field theories, the Formulae (20)–(23) have to be supplemented by the respective transformation formulae for the fields in order to define the transformations completely. For demonstration purposes consider the Lagrangian (16) for inviscid barotropic flows in the absence of external forces, *U* = 0:

$$\ell\left(\varrho, a, \dot{\varrho}, \nabla\varrho, \dot{\beta}, \nabla\beta\right) = -\varrho \left[\dot{\varrho} + a\dot{\beta} + \frac{1}{2}\left(\nabla\varrho + a\nabla\beta\right)^2 + e\left(\varrho\right)\right],\tag{24}$$

where the dot indicates the partial time derivative *∂*/*∂t*. Obviously, the Lagrangian (24) is invariant with respect to space and time translations in Equations (20), (21) and rigid rotations (22) if the four fields, , *ϕ*, *α*, *β* are assumed to be likewise invariant. In contrast, invariance with respect to Galilei boosts requires the first Clebsch variable *ϕ* to be transformed according to [26]:

$$
\varphi \to \varphi' = \varphi - \vec{u}\_0 \cdot \vec{x} + \frac{\vec{u}\_0^2}{2}t,\tag{25}
$$

implying, according to Equation (5), that consequently *u* → *u* = *u* − *u*<sup>0</sup> for the velocity field and compensating the non-invariance *∂*/*∂t* → *∂*/*∂t* + *u*<sup>0</sup> · ∇ of the partial time derivatives occurring in the Lagrangian. The fact that the Galilei boost becomes manifest as a combination of the geometric transformation in Equation (23) with the gauge transformation in Equation (25) is an unfavourable feature since, for each physical system depending on an arbitrary set of fields, the transformation formulae for the fields involved have to be defined individually. However, there is a way to eliminate this disadvantage: by means of the substitution:

$$
\varphi = \Phi + \mathcal{J},
\tag{26}
$$

with the *generating field*:

$$
\zeta := \frac{\vec{x}^2}{2t}.\tag{27}
$$

The Lagrangian (24) can be re-written in the alternative form:

$$\tilde{\ell}\left(\varrho, a, \stackrel{\circ}{\Phi}, \nabla\Phi, \stackrel{\circ}{\beta}, \nabla\beta\right) = -\varrho\left[\stackrel{\circ}{\Phi} + a\stackrel{\circ}{\beta} + \frac{1}{2}\left(\nabla\varphi + a\nabla\beta\right)^2 + e\left(\varrho\right)\right],\tag{28}$$

where the ring symbol ◦ indicates the *dual time derivative* [56]:

$$\stackrel{\circ}{\Phi} = \left\{ \frac{\partial}{\partial t} + \nabla \zeta \cdot \nabla \right\} \Phi \,, \tag{29}$$

which in contrast to the conventional time derivative is invariant with respect to Galilei boosts. As a consequence, the Lagrangian in its alternative form in Equation (28), subsequently called the *dual representation*, proves to be invariant with respect to Galilei boosts if all fields including Φ are assumed to be likewise invariant. Thus, the essence of the dual representation in Equation(28) is that Galilei boosts become manifest as pure geometrical transformations without the need to combine them with a re-gauging of potentials. Vice versa, for the dual Lagrangian (28), time and space translations become manifest as mixed transformations consisting of a geometric part in combination with a re-gauging of the potential Φ [56]. Subsequently, a Lagrangian is termed *strictly invariant* if it is invariant with respect to a purely geometric transformation without the need of re-gauging one of the potential fields. In this sense, Equation (24) is strictly invariant with respect to space and time translations while Equation (28) is strictly invariant with respect to Galilei boosts but not vice versa.

This poses the question of whether the coexistence of two different representations of a given Lagrangian, one being strictly invariant with respect to space and time translations and the other one being strictly invariant with respect to Galilei boosts, is an individual feature of the theory of inviscid barotropic flows or if other examples exist. Scholle [56] undertook a rigorous analysis to investigate this question, arriving at the conclusion that this coexistence is a *universal property of every continuum* being ruled by self-adjoint equations of motions invariant with respect to the full Galilean group: if the state of a continuous system is defined by *N* independent fields *ψ<sup>i</sup>* , *i* = 1 ··· , *N* and its evolution determined by a variational principle:

$$\delta \int\_{t\_1}^{t\_2} \iiint\limits\_V \ell \left(\psi\_{i\prime} \psi\_{i\prime} \nabla \psi\_i\right) \mathrm{d}V \mathrm{d}t = 0 \,,\tag{30}$$

based on free and independent variation of the fields *ψ<sup>i</sup>* and their first order spatial and temporal derivatives for fixed values at the initial and final time, *t*1,2, then a variable transformation

$$
\Psi\_i = \mathbb{K}\_i \left( \mathbb{Y}\_{j\prime} \mathbb{Z}\_{\prime\prime} \nabla \mathbb{Q} \right) \,, \tag{31}
$$

exists with the generating field *ζ* defined via Equation (27), fulfilling the identity:

$$\ell\left(\mathbb{K}\_{i}\left(\mathbb{Y}\_{j},\mathbb{Y}\_{\prime}\nabla\mathbb{S}\right),\dot{\mathbb{K}}\_{i}\left(\mathbb{Y}\_{j},\mathbb{Y}\_{\prime}\nabla\mathbb{S}\right),\nabla\mathbb{K}\_{i}\left(\mathbb{Y}\_{j},\mathbb{Y}\_{\prime}\nabla\mathbb{S}\right)\right) = \ell\left(\mathbb{Y}\_{i},\overset{\circ}{\mathbb{Y}}\_{i},\nabla\mathbb{Y}\_{i} + \frac{1}{t}\overrightarrow{\mathbb{K}}\_{i}\left(\mathbb{Y}\_{j}\right)\right),\tag{32}$$

$$\mathcal{K}\_i\left(\Psi\_j\right) := \lim\_{\substack{\zeta, \nabla \zeta \to 0 \ \overline{\partial}\left(\nabla \zeta\right) \text{ } }} \frac{\partial \mathcal{K}\_i}{\partial \left(\nabla \zeta\right) \text{ } }\, \, \, \, \tag{33}$$

giving rise to the dual representation of the Lagrangian on the right hand of Equation (32) depending on the dual time derivatives ◦ <sup>Ψ</sup>*i*<sup>=</sup> <sup>Ψ</sup>˙ *<sup>i</sup>* <sup>+</sup> <sup>∇</sup>*<sup>ζ</sup>* · ∇Ψ*<sup>i</sup>* of the transformed fields.

Since the conventional representation (*ψi*, *ψ*˙ *<sup>i</sup>*, ∇*ψi*) is obviously strictly invariant with respect to space and time translations while the dual representation Ψ*i*, ◦ <sup>Ψ</sup>*i*, <sup>∇</sup>Ψ*<sup>i</sup>* <sup>+</sup> *<sup>K</sup> <sup>i</sup>* Ψ*j* /*t* is strictly invariant with respect to Galilei boosts, simultaneous invariance with respect to translations and Galilei boosts is granted by Equation (32), which consequentially can be understood as a collective symmetry criterion for the Galilean group and has to be fulfilled by any Lagrangian related to a physically closed continuous system. One implication is that the gauge transformation:

$$
\psi\_i \to \psi\_i' = K\_i \left( \psi\_{\bar{\jmath}'}, \varepsilon, 0 \right) \tag{34}
$$

is likewise defined as being a symmetry transformation of the Lagrangian. This *induced gauge transformation* is remarkable since it is an indirect consequence of the Galilean invariance.

The above-mentioned general properties of Lagrangians in continuum theory entail additional general implications for the physical balances resulting from the variational Principle (30) by utilising Noether's theorem [57–59], the essence of which is that each Lie symmetry of the Lagrangian gives rise to a physical balance equation and to a canonical definition of the density and flux density involved. In the present context, the balances for mass and momentum,

$$\frac{\partial \varrho}{\partial t} + \nabla \cdot (\varrho \vec{u}) = 0 \,, \tag{35}$$

$$
\frac{
\partial \vec{p}
}{
\partial t
} + \nabla \cdot \underline{\Pi} = \vec{0} \; \text{,}
\tag{36}
$$

are related (via Noether's theorem) to the induced gauge Transformation (34) and (from the space translation in Equation (21)) with the mass density , the mass flux density*j* = *u*, the momentum density *p* and the momentum flux density Π which are given by [56]:

$$\varrho = -\frac{\Im \ell}{\partial \dot{\varphi}\_i} K\_{0i} \left( \psi\_{\dot{\jmath}} \right) \ . \tag{37}$$

$$
\vec{j} = \varrho \vec{u} = -\frac{\partial \ell}{\partial \nabla \psi\_i} \mathcal{K}\_{0i} \left( \psi\_j \right) \; , \tag{38}
$$

$$
\vec{p} = -\frac{\partial \ell}{\partial \dot{\psi}\_i} \nabla \psi\_i \,\,\,\,\,\tag{39}
$$

$$
\underline{\Pi} = \ell \underline{1} - \frac{\partial \ell}{\partial \nabla \psi\_i} \otimes \nabla \psi\_i \, \, \, \, \tag{40}
$$

with the infinitesimal Generator *K*0*<sup>i</sup> ψj* defined as:

$$\mathcal{K}\_{0i} \left( \Psi\_j \right) := \lim\_{\substack{\mathcal{I}\_\star, \nabla \mathcal{I}\_\star \to 0 \ \partial \mathcal{I}\_\star}} \frac{\partial \mathcal{K}\_i}{\partial \mathcal{J}\_\star} \,. \tag{41}$$

In Equations (37)–(40) and subsequently, the following *Einstein notation* is utilised: an index variable appearing twice in a single term implies summation of that term over all values of the index.

Note that in contrast to classical continuum mechanics, the mass flux density and the momentum density need not to be equal. The difference between both, *<sup>p</sup>* <sup>∗</sup> :<sup>=</sup> *<sup>p</sup>* <sup>−</sup>*<sup>j</sup>* <sup>=</sup> *<sup>p</sup>* <sup>−</sup> *u*, is called the *quasi-momentum density* and can be interpreted as contributions to the momentum due to non-material degrees of freedom, for example, electromagnetic fields, thermal fields and also due to phenomena beyond the scope of continuum theories on a molecular scale, for example, Brownian motion. After analysing the Noether balance resulting from Galilei-boosts, a constitutive relation between momentum and mass flux can be identified [56] implying the identity:

$$\vec{p}^{\*} = -\frac{\partial}{\partial t} \left[ \frac{\partial \ell}{\partial \psi\_{i}} \vec{\mathcal{K}}\_{i} \left( \psi\_{j} \right) \right] - \nabla \cdot \left[ \frac{\partial \ell}{\partial \nabla \psi\_{i}} \otimes \vec{\mathcal{K}}\_{i} \left( \psi\_{j} \right) \right] \tag{42}$$

for the quasi-momentum density.

In classical continuum mechanics, mass flux density and momentum density are assumed to be equal, implying *p* <sup>∗</sup> = 0, which according to Equations (42) and (33) requires the variable Transformation (31) to be independent of ∇*ζ*, leading to a drastic mathematical simplification of Criterion (32) and allowing derivation of the following universal scheme for Lagrangians [56]: without loss of generality the conventional representation of the Lagrangian can be written in terms of *N* fields, (*ψi*) = (*ϕ*, *ϑ*1, ··· , *ϑN*−1) as:

$$\ell\left(\dot{\mathfrak{o}}, \nabla\,\mathfrak{o}, \mathfrak{o}^{\dot{\mathfrak{o}}}, \dot{\mathfrak{o}}^{\dot{\mathfrak{o}}}, \nabla\,\mathfrak{o}^{\dot{\mathfrak{o}}}\right) = L\left(\omega, \mathfrak{o}\_{\dot{\mathfrak{o}}}, \dot{\mathfrak{o}}\_{\dot{\mathfrak{o}}}, \nabla\mathfrak{o}\_{\dot{\mathfrak{o}}}\right),\tag{43}$$

with the abbreviations:

$$
\omega := \phi + \frac{1}{2} \left( \nabla \phi \right)^2 \,, \tag{44}
$$

$$\dot{\hat{\theta}}\_{\dot{j}}^{\circ} := \dot{\theta}\_{\dot{j}} + \nabla \boldsymbol{\varphi} \cdot \nabla \theta\_{\dot{j}}, \qquad \dot{j} = 1, \cdots, N - 1 \; , \tag{45}$$

fulfilling Criterion (32) for the following form of the variable Transformation (31):

$$
\varphi = \Phi + \zeta \,, \tag{46}
$$

$$
\theta\_{\rangle} = \Theta\_{\rangle} \,, \qquad j = 1, \cdots, N - 1 \; . \tag{47}
$$

Consequently, the Lagrangian (24) for inviscid barotropic flows corresponds to the universal Scheme (43) with *L* = − *ω* + *α β* +<sup>1</sup> <sup>2</sup> (*α*∇*β*) <sup>2</sup> + *e*() but also any other Lagrangian of a continuous system with Galilean invariance and equality of momentum density and mass flux density. Moreover, in [56], it is certified that the velocity field fulfils all transformation rules with respect to the Galilean Group in any case.

For the universal Scheme (43), the mass density (Equation ((37)) simplifies to

 <sup>=</sup> <sup>−</sup> *<sup>∂</sup><sup>L</sup> ∂ω* , (48)

while the momentum density, which equals the mass flux density, according to Equation (39), becomes:

$$\delta\vec{\mu} = \vec{p} = \underbrace{-\frac{\partial L}{\partial \omega}}\_{\vec{\Theta}} \nabla \varphi - \frac{\partial L}{\partial \left. \vec{\theta} \right|\_{i}} \nabla \theta\_{i} \, \, \, \tag{49}$$

Therefore, the velocity field takes the form:

*u* = ∇*ϕ* + *γi*∇*ϑ<sup>i</sup>* , (50)

of a generalised Clebsch representation with:

$$\gamma\_i = \gamma\_i \left( \omega, \theta\_{\dot{\jmath}}, \stackrel{\circ}{\theta}\_{\dot{\jmath}}, \nabla \theta\_{\dot{\jmath}} \right) := -\frac{1}{\varrho} \frac{\partial L}{\partial \stackrel{\circ}{\theta}\_{\dot{\jmath}}}.\tag{51}$$

The surprising result from the above investigations is that the Clebsch representation of the velocity field is an inevitable consequence of the Galilean invariance. More precisely, the essence of the above can be formulated in the following theorem:

**Theorem 1.** *If a continuum, no matter whether it is a solid or a fluid, fulfils the following three requirements:*


*then the velocity field takes, via Noether's theorem, the analytic form (Equation* (50)*) of a (generalised) Clebsch representation.*

In terms of the above theorem, the Clebsch representation is substantiated from fundamental symmetries in physics.

#### *2.3. An Extended Clebsch Transformation for Viscous Flow*

Consider now the NS equations together with the continuity equation:

$$
\frac{\text{D}\vec{u}}{\text{D}t} - \nu \Delta \vec{u} + \nabla \left[\frac{p}{\varrho} + \mathcal{U}\right] \quad = \quad \vec{0}, \tag{52}
$$

$$
\nabla \cdot \vec{u} \, \, \, \, = \, \, 0 \,, \, \, \tag{53}
$$

assuming incompressible flow according to Equation (53), commensurate with the classical theory of viscous flow [19], such that *ν* denotes the kinematic viscosity. As demonstrated in Section 2.1 for inviscid flow, the key feature of the Clebsch transformation is that it enables the Euler equations to be written in the form of Equation (10), being a natural outcome of the mathematical structure of inviscid

*Water* **2020**, *12*, 1241

flow theory. For the viscous case, it is less obvious how the transformation applies, since the specific viscous force −*ν*Δ*u* in the NS equations reads, when written in terms of the Clebsch variables, as:

$$-\nu\Delta\vec{u} = 2\nu\nabla \times \vec{w} - 2\nu\nabla\underbrace{(\nabla \cdot \vec{u})}\_{0} = \nu\Delta\p\nabla\vec{u} - \nu\Delta\vec{a}\nabla\beta - \nu\left(\nabla\vec{a}\cdot\nabla\right)\nabla\beta + \nu\left(\nabla\beta\cdot\nabla\right)\nabla\vec{a} \,, \tag{54}$$

which does not correspond directly to the mathematical Scheme (10). Despite this, it is shown in [60] that it is possible to decompose the viscous term according to Equation (10) by introducing an additional field *ξ*. In addition and as demonstrated subsequently, this procedure applies to any vector field, not only to the specific viscous friction force (Equation (54)).

**Theorem 2.** *Leta be an arbitrary vector field, u the velocity field given in Clebsch representation according to Equation* (5) *based on the Clebsch variables ϕ*, *α*, *β and ω the vorticity given according to Equation* (14)*. Then, a decomposition ofa according to Scheme* (10) *is always possible, i.e., three fields ξ, λ and μ exist fulfilling:*

$$
\overrightarrow{a} = \nabla \xi + \lambda \nabla \pi + \mu \nabla \beta \,\tag{55}
$$

*for prescribed fields α*, *β where the auxiliary field ξ results as the solution of the conditional equation:*

$$
\vec{\omega} \cdot \nabla \vec{\xi} = \vec{\omega} \cdot \vec{\pi} \tag{56}
$$

*and the two coefficients λ*, *μ explicitly as:*

$$
\lambda = \frac{\vec{\omega} \times [\vec{a} - \nabla \vec{\xi}]}{2\vec{\omega}^2} \cdot \nabla \beta \,\tag{57}
$$

$$
\mu = -\frac{\vec{\omega} \times \left[\vec{\mu} - \nabla \vec{\xi}\right]}{2\vec{\omega}^2} \cdot \nabla \vec{\alpha} \,. \tag{58}
$$

**Proof.** The conditional Equation (56) is obtained by taking the inner product of Equation (55) with 2*ω* = ∇*α* × ∇*β*. Once having solved Equation (56), the identity:

$$
\vec{\omega} \times \left( \vec{\omega} \times \left[ \vec{a} - \nabla \vec{\xi} \right] \right) = \vec{\omega} \underbrace{\left( \vec{\omega} \cdot \left[ \vec{a} - \nabla \vec{\xi} \right] \right)}\_{0} - \left[ \vec{a} - \nabla \vec{\xi} \right] \vec{\omega}^{2} = -\vec{\omega}^{2} \left[ \vec{a} - \nabla \vec{\xi} \right] \tag{59}
$$

follows and therefore:

$$\begin{split} \vec{\omega} - \nabla \xi &= \frac{(\vec{\omega} \times [\vec{\omega} - \nabla \xi]) \times \vec{\omega}}{\vec{\omega}^2} = \underbrace{\frac{(\vec{\omega} \times [\vec{\pi} - \nabla \xi]) \times (\nabla a \times \nabla \beta)}{2\vec{\omega}^2}}\_{=:\lambda} \\ &= \underbrace{\frac{(\vec{\omega} \times [\vec{\pi} - \nabla \xi]) \cdot \nabla \beta}{2\vec{\omega}^2}}\_{=:\lambda} \nabla a \underbrace{-\frac{(\vec{\omega} \times [\vec{\pi} - \nabla \xi]) \cdot \nabla a}{2\vec{\omega}^2}}\_{=:\mu} \nabla \beta \end{split} \tag{60}$$

as the desired decomposition of *a* − ∇*ξ* as a linear combination of ∇*α* and ∇*β* according to Equation (55).

Like the Clebsch variables Φ, *α*, *β*, the auxiliary field *ξ* is not uniquely given, since any particular solution *ξ*<sup>p</sup> of the inhomogeneous linear first order PDE (56) can be superposed with any solution *ξ*<sup>h</sup> of the respective homogeneous PDE *ω* · ∇*ξ*<sup>h</sup> = 0. Since three independent solutions are given by *α*, *β* and *t*, the mathematical theory of linear first order PDE implies *ξ*<sup>h</sup> = *F*(*α*, *β*, *t*) for an arbitrary function *F*. As a consequence:

$$
\mathfrak{J} \longrightarrow \mathfrak{J}' = \mathfrak{J} + F(\mathfrak{a}, \mathfrak{b}, t),
\tag{61}
$$

is a gauge transformation for the auxiliary field, which is used subsequently to provide a favourable form of the resulting equations.

*Water* **2020**, *12*, 1241

Note that the decomposition in Equation (60) is applicable to an arbitrary vector field *a*, including for instance non-conservative forces. In the case of viscous flow, the specific viscous force in Equation (54) in the NS Equations (52) can be partitioned into two parts [60], namely *<sup>a</sup>* = *<sup>ν</sup>* (∇*<sup>β</sup>* · ∇) ∇*<sup>α</sup>* − *<sup>ν</sup>* (∇*<sup>α</sup>* · ∇) ∇*<sup>β</sup>* on the one hand and −*ν*∇2*<sup>u</sup>* −*<sup>a</sup>* = *<sup>ν</sup>*Δ*β*∇*<sup>α</sup>* − *<sup>ν</sup>*Δ*α*∇*<sup>β</sup>* on the other, which makes sense at first glance since the latter contributions already have the desired form of Equation (10). However, a major disadvantage of this partition is that the selected field *a* is not gauge invariant with respect to the transformation shown in Equation (6). Although not representing a problem mathematically, this feature makes the field equations resulting from the generalised Clebsch transformation less transparent from a physical viewpoint.

In contrast to the above, the full specific viscous force in Equation (54) is considered in the following, i.e.,

$$
\vec{a} = 2\nu \nabla \times \vec{\omega} \, , \tag{62}
$$

leading to a physically more transparent representation of the field equations, as provided in [60]. This procedure allows the NS Equations (52) to be written as:

$$\vec{\nabla} = \nabla \left[ \frac{\partial \underline{\rho}}{\partial t} + \kappa \frac{\partial \underline{\beta}}{\partial t} + \frac{\vec{u}^2}{2} + \frac{p}{\underline{\rho}} + \mathcal{U} + \vec{\xi} \right] + \left[ \frac{\mathcal{D}a}{\mathcal{D}t} + \mu \right] \nabla \beta + \left[ -\frac{\mathcal{D}\beta}{\mathcal{D}t} + \lambda \right] \nabla a \,, \tag{63}$$

with *ξ*, *λ*, *μ* given according to Theorem 2. As in Section 2.1, by proper gauging of the potentials, the bracketed terms in Equation (63) vanish separately, resulting in the following set of PDEs

$$
\frac{\partial \mathcal{Q}}{\partial t} + a \frac{\partial \mathcal{B}}{\partial t} + \frac{i\vec{l}^2}{2} + \frac{p}{\varrho} + l\mathcal{I} + \xi = 0\tag{64}
$$

$$\frac{\text{D}\mathfrak{a}}{\text{D}t} - \frac{\vec{\omega} \times \left[2\nu \nabla \times \vec{\omega} - \nabla \xi\right]}{2\vec{\omega}^2} \cdot \nabla a = 0 \tag{65}$$

$$\frac{\text{D}\beta}{\text{D}t} - \frac{\vec{\omega} \times [2\nu \nabla \times \vec{\omega} - \nabla \xi]}{2\vec{\omega}^2} \cdot \nabla \beta = 0 \tag{66}$$

supplemented by the conditional equation:

$$
\vec{\omega} \cdot \nabla \xi = 2\nu \vec{\omega} \cdot (\nabla \times \vec{\omega}) \tag{67}
$$

for the auxiliary field and the continuity Equation (53). The latter, in terms of Clebsch variables [43], becomes:

$$
\nabla^2 \varrho + \mathfrak{a} \nabla^2 \beta + 2 \nabla \mathfrak{a} \cdot \nabla \beta = 0 \,. \tag{68}
$$

Physically, Equation (64) can be interpreted as a generalised Bernoulli's equation and Equations (65) and (66) as generalised transport equations covering the entire vortex dynamics of the flow.

#### *2.4. Axisymmetric Stagnation Flow*

The problem of an axisymmetric stagnation flow against a solid wall, see Figure 1, is considered one of prototypical character in fluid mechanics, since the underlying features of Navier–Stokes theory are exposed, especially the formation of a boundary layer at the solid wall *z* = 0. It is one of the few boundary layer problems that allow for an analytical treatment of the full NS equations, without the necessity of neglecting higher order terms. In the inviscid case, the velocity field, written in cylindrical coordinates (*r*, *φ*, *z*), is given by Mayes et al. [61]:

$$
\vec{u}\_{\text{invis}} = A r \vec{e}\_{\text{I}} - 2A z \vec{e}\_{z} \,. \tag{69}
$$

Although Equation (69) is a solution of the NS equations, it does not match to the no-slip condition *er* ·*u* = 0 at the wall and therefore represents a good approximation for the far field only, i.e., the field in the region *<sup>z</sup>* <sup>√</sup>*ν*/*<sup>A</sup>* far beyond the boundary layer; whereas, in the vicinity of the wall, a boundary layer becomes manifest in which the velocity is assumed to take the slightly different form:

$$\vec{u} = rf'(z)\vec{e}\_r - 2f(z)\vec{e}\_z\,. \tag{70}$$

containing the function *f*(*z*) that has to be determined. Note that the continuity Equation (53) is fulfilled identically by Equation (70). The associated boundary conditions are (i) the no-slip/no-penetration condition *<sup>u</sup>* <sup>=</sup>0 at *<sup>z</sup>* <sup>=</sup> 0 and (ii) the matching condition *<sup>u</sup>* <sup>→</sup> *u*invis for *<sup>z</sup>* <sup>→</sup> <sup>∞</sup>, leading to:

$$f(0) = 0, \qquad f'(0) = 0, \qquad \lim\_{z \to \infty} f'(z) = A. \tag{71}$$

Now the extended Clebsch transformation developed in Section 2.3 is applied: writing the vorticity and the specific viscous force as:

$$
\vec{\omega} = \frac{r}{2} f''(z) \,\vec{e}\_{\phi} \,\prime \tag{72}
$$

$$2\nu \nabla \times \vec{\omega} = \nu \left[ 2f^{\prime\prime}(z) \, \vec{e}\_z - rf^{\prime\prime}(z) \, \vec{e}\_r \right] \, \,, \tag{73}$$

the conditional equation simplifies to:

$$\frac{\partial \xi}{\partial \phi} = 0,\tag{74}$$

implying *ξ* = *ξ* (*r*, *z*, *t*) as a general solution. Since a steady axisymmetric flow is considered, it is subsequently assumed independent of time and azimuthal angle for all potentials *ξ*, *α*, *β* and *ϕ*. The remaining equations of interest are the generalised transport Equations (65) and (66) which, after minor mathematical manipulation, take the form:

$$\left[r^2 f'\left(z\right) f''\left(z\right) - 2\nu f''\left(z\right) + \frac{\partial \xi}{\partial z}\right] \frac{\partial a}{\partial r} - \left[2rf\left(z\right) f''\left(z\right) + \nu r f'''\left(z\right) + \frac{\partial \xi}{\partial r}\right] \frac{\partial a}{\partial z} = 0,\tag{75}$$

$$\left[r^2 f'(z) f'''(z) - 2\nu f'''(z) + \frac{\partial \xi}{\partial z}\right] \frac{\partial \beta}{\partial r} - \left[2rf(z) f'''(z) + \nu rf''''(z) + \frac{\partial \xi}{\partial r}\right] \frac{\partial \beta}{\partial z} = 0. \tag{76}$$

The generalised Bernoulli Equation (64) and the continuity Equation (68) are decoupled from Equations (75) and (76); they provide the third Clebsch variable *ϕ* and the pressure *p*, both of which are not required here.

Three unknown fields *ξ*, *α*, *β* are involved in the two decisive PDEs (75) and (76) due to the freedoms given by the gauge symmetry with respect to the transformation in Equation (6). This provides the opportunity to choose the potentials in a beneficial way: for instance, by considering that the vorticity in Equation (72) can be written as 2*ω* = *r f* (*z*)*e<sup>φ</sup>* = *r f* (*z*)*ez* ×*er* = ∇ *<sup>f</sup>* (*z*) × ∇*r*2/2 and comparison with Equation (14) motivates the choice:

$$
\mathfrak{a} = f'(z) \tag{77}
$$

$$
\beta = \frac{r^2}{2} \tag{78}
$$

for the two Clebsch variables. By inserting these into Equations (75) and (76), the following simplified equations:

$$\frac{\partial \overline{\xi}}{\partial r} = -2rf(z)f''(z) - \nu rf''''(z) \ , \tag{79}$$

$$\frac{\partial \bar{\xi}}{\partial z} = 2\nu f''\left(z\right) - r^2 f'\left(z\right) f''\left(z\right) \, , \tag{80}$$

result for the auxiliary field which are both integrable, leading to:

$$\xi = -r^2 f(z) f''''(z) - \frac{\nu}{2} r^2 f''''(z) + F\_1(z) \, . \tag{81}$$

$$
\xi = 2\nu f'\left(z\right) - \frac{r^2}{2} f'\left(z\right)^2 + F\_2(r) \,. \tag{82}
$$

containing the two integration functions *F*1(*z*) and *F*2(*r*). By subtracting Equation (81) from (82), the auxiliary field is eliminated, giving:

$$\frac{\sigma^2}{2} \left[ \nu f^{\prime\prime\prime}(z) + 2f(z)f^{\prime\prime}(z) - f^{\prime}(z)^2 \right] + F\_2(r) + 2\nu f^{\prime}(z) - F\_1(z) = 0. \tag{83}$$

Next, taking the limit *r* → 0 of the above equation leads to *F*1(*z*) = *F*2(0) + 2*ν f* (*z*); alternatively, the limit *z* → ∞, in tandem with the matching condition *f* (*z*) → *A* according to Equation (71), gives *<sup>F</sup>*2(*r*) − *<sup>F</sup>*2(0) = *<sup>A</sup>*2*r*2/2. On re-inserting these explicit forms of the integration functions into Equation (83), the following ODE for the unknown function *f*(*z*) is obtained:

$$\nu f''''(z) + 2f(z)f''(z) - f'(z)^2 + A^2 = 0. \tag{84}$$

Finally, the substitution *<sup>f</sup>*(*z*) = <sup>√</sup>*ν<sup>A</sup>* ¯ *<sup>f</sup>*(*z*¯) with *<sup>z</sup>*¯ <sup>=</sup> <sup>√</sup>*A*/*ν<sup>z</sup>* transforms the latter into the more convenient form:

$$
\bar{f}'''(\bar{z}) + 2\bar{f}(\bar{z})\bar{f}''(\bar{z}) - \bar{f}'(\bar{z})^2 + 1 = 0,\tag{85}
$$

well-known in the literature as the *Falkner-Skan equation* [61].

**Figure 1.** Schematic of an axisymmetric stagnation flow in the vicinity of a solid wall.

#### **3. Complex Variable and Tensor Potential Approach**

#### *3.1. The Classical Complex Variable Method*

Consider the case of an arbitrary homogeneous continuum, with plane stress given by the tensor:

$$
\underline{T} = \begin{pmatrix} \sigma\_{\rm x} & \tau\_{\rm xy} \\ \tau\_{\rm xy} & \sigma\_{\rm y} \end{pmatrix} \tag{86}
$$

under a load provided by a conservative force with specific potential energy *U*. Assuming a steady state, the equilibrium condition:

$$
\nabla \cdot \underline{T} + \varrho \nabla \iota \underline{I} = \vec{0} \tag{87}
$$

has to be fulfilled. By defining the complex variable:

$$\mathbf{y} := \mathbf{x} + \mathbf{i}y \,, \tag{88}$$

*Water* **2020**, *12*, 1241

the three fields *σx*, *σ<sup>y</sup>* and *τxy* can be considered as functions of *ξ* and its complex conjugate, *ξ*, and the equilibrium Condition (87) reads:

$$\begin{split} \vec{0} &= \left( \frac{\partial}{\partial \xi} + \frac{\partial}{\partial \xi}, \mathbf{i}\frac{\partial}{\partial \xi} - \mathbf{i}\frac{\partial}{\partial \xi} \right) \left( \begin{array}{c} \sigma\_{\text{x}} + \varrho \boldsymbol{l} \boldsymbol{I} \\ \tau\_{\text{xy}} \end{array} \begin{array}{c} \tau\_{\text{xy}} \\ \sigma\_{\text{y}} + \varrho \boldsymbol{l} \boldsymbol{I} \end{array} \right) \\ &= \frac{1}{2} \left( \begin{array}{c} \frac{\partial}{\partial \xi} \left[ \sigma\_{\text{x}} + \sigma\_{\text{y}} + 2\varrho \boldsymbol{l} \boldsymbol{I} \right] + \frac{\partial}{\partial \xi} \left[ \sigma\_{\text{x}} - \sigma\_{\text{y}} + \mathrm{i}\tau\_{\text{xy}} \right] + \mathrm{c.c.} \\ -\mathrm{i}\frac{\partial}{\partial \xi} \left[ \sigma\_{\text{x}} + \sigma\_{\text{y}} + 2\varrho \boldsymbol{l} \boldsymbol{I} \right] - \mathrm{i}\frac{\partial}{\partial \xi} \left[ \sigma\_{\text{x}} - \sigma\_{\text{y}} + \mathrm{i}\tau\_{\text{xy}} \right] + \mathrm{c.c.} \end{array} \right) \end{split} \tag{89}$$

where c.c. denotes the conjugate complex of the preceding expression. It is convenient to introduce the hydrostatic pressure *p* := (*σ<sup>x</sup>* + *σy*)/2 and the complex stress *σ* := (*σ<sup>x</sup>* − *σy*)/2 + i*τxy*, allowing Equation (89) to be written as *<sup>∂</sup>* [*<sup>p</sup>* + *U*] /*∂ξ* + *∂σ*/*∂ξ* = 0 and *∂* [*p* + *U*] /*∂ξ* + *∂σ*/*∂ξ* = 0 and therefore:

$$
\frac{\partial}{\partial \overline{\xi}} \left[ \frac{p}{\varrho} + \mathcal{U} \right] + \frac{1}{\varrho} \frac{\partial \underline{\varphi}}{\partial \underline{\xi}} = 0,\tag{90}
$$

as the complex form of the equilibrium Condition (87). The equilibrium Condition (90) is fulfilled identically by introducing a real-valued potential field Φ according to:

$$\frac{\mathcal{\underline{\mathcal{F}}}}{\mathcal{\underline{\mathcal{G}}}} = -4 \frac{\partial^2 \Phi}{\partial \underline{\mathcal{G}}^2} \, , \tag{91}$$

$$\frac{p}{\varrho} + \mathcal{U} = 4 \frac{\partial^2 \Phi}{\partial \overline{\xi^2} \partial \overline{\xi}}. \tag{92}$$

The potential Φ is the well known *Airy stress function* from the theory of linear elasticity [32,62]. Since, however, no assumptions regarding the constitutive equations of the respective continuum are required for the above derivation of the complex equilibrium Condition (87) and the use of Airy's stress function according to Equations (91) and (92), this approach applies to any continuum and thus also to Stokes flow. By assuming a steady flow and neglecting the nonlinear terms for the case of very small Reynolds numbers, the NS Equations (52) simplify to the Stokes equation; being of the general form in Equation (87) with stress tensor *Tij* = −*pδij* + *ν ∂iuj* + *∂jui* , leading to a complex stress field: *<sup>σ</sup>*

$$\frac{\mathcal{Q}}{\mathcal{\varrho}} = -2\nu \frac{\partial \mu}{\partial \mathcal{\xi}'} \tag{93}$$

where *u* denotes the complex velocity:

$$
\mu = \mu\_{\mathbf{x}} + \mathbf{i}\mu\_{\mathbf{y}}.\tag{94}
$$

On introducing a stream function Ψ, satisfying:

$$
\mu = -2\mathrm{i}\frac{\partial\Psi}{\partial\tilde{\xi}}\,'\,\tag{95}
$$

which fulfils the continuity Equation (53) identically [37,38], the constitutive Equation (93), written in terms of the stream function, becomes:

$$\frac{\varrho}{\varrho} = 4\text{i}\nu \frac{\partial^2 \Psi}{\partial \tilde{\xi}^2}. \tag{96}$$

By inserting Equation (96) into (91), the fully integrable equation:

$$\frac{\partial^2}{\partial \overline{\zeta}^2} \left[ \Phi + \mathrm{i} \nu \Psi \right] = 0,\tag{97}$$

results; the double integration of which gives the general solution:

$$
\Phi + \mathrm{i}\nu\Psi = \mathrm{g}\_0\left(\xi\right) + \xi\mathrm{g}\_1\left(\xi\right),
\tag{98}
$$

containing integration functions *g*<sup>0</sup> (*ξ*) and *g*<sup>1</sup> (*ξ*), known as *Goursat functions* [62]. Note that, such an approach leads to the well-known Sherman–Lauricella equations [63–65] and to a reproduction of the Muskhelishvili–Kolosov formula [32,33].

The complex variable method has been applied successfully to various Stokes' flow problems [35,66–69], usually by adopting a conformal mapping on the Goursat functions known for a simple flow domain in order to obtain the respective Goursat functions for a non-trivial domain. Being holomorphic, the Goursat functions can be recovered alternatively from their boundary values determined by the boundary conditions as utilised in [70–72] as an investigative tool for studying the internal flow structure of film flows over corrugated walls.

#### *3.2. Integration of the Full 2D Navier–Stokes Equations*

Making use again of the complex variable transformation in Equations (88) and (94), the NS Equations (52) and continuity Equation (53) can be reformulated as:

$$2\frac{\partial\mu}{\partial t} + 2\frac{\partial}{\partial\xi}\left[\frac{\bar{u}\mu}{2} + \frac{p}{\varrho} + \mathcal{U}\right] + 2\frac{\partial}{\partial\xi}\left(\frac{\mu^2}{2}\right) = 4\nu\frac{\partial^2\mu}{\partial\xi^2\partial\xi'}\tag{99}$$

$$\Re\left(\frac{\partial u}{\partial \xi}\right) = 0,\tag{100}$$

in terms of the complex coordinate *ξ*, the complex velocity field *u* and their complex conjugates, with denoting the real part of the subsequent complex expression. It is obvious that by introducing a stream function Ψ, according to Equation (95), the continuity Equation (100) is fulfilled identically. Accordingly, the complex NS Equation (99) can be written as:

$$\frac{\partial}{\partial \overline{\xi}} \left[ \frac{\partial u}{2} + \frac{p}{\varrho} + \mathcal{U} - \mathbf{i} \frac{\partial \Psi}{\partial t} \right] + \frac{\partial}{\partial \overline{\xi}} \left( \frac{u^2}{2} \right) = 2\nu \frac{\partial^2 u}{\partial \overline{\xi} \partial \overline{\xi}}.\tag{101}$$

By introducing a new complex potential *M* according to

$$\frac{\partial \mathcal{U}}{2} + \frac{p}{\varrho} + \mathcal{U} - \mathrm{i}\frac{\partial \Psi}{\partial t} = 2\frac{\partial \mathcal{M}}{\partial \underline{\chi}}\,'\,. \tag{102}$$

an integrable form of Equation (99),

$$2\frac{\partial}{\partial\xi}\left[2\frac{\partial M}{\partial\xi^T} + \frac{u^2}{2} - 2\nu\frac{\partial u}{\partial\xi^T}\right] = 0 \,,\tag{103}$$

is obtained which, after integration with respect to *ξ*, gives:

$$2\frac{\partial M}{\partial \overline{\xi}^{\overline{\eta}}} + \frac{\mu^2}{2} - 2\nu \frac{\partial \mu}{\partial \overline{\xi}^{\overline{\eta}}} = f(\overline{\xi}),\tag{104}$$

having integration function *f*(*ξ*) on the right hand side. The latter can be conveniently set to zero by re-gauging the potential *M* as follows:

$$M \longrightarrow M + \frac{1}{2}F(\overline{\xi}),$$

with *F* (*ξ*) = *f*(*ξ*), since according to Equation (102) any complex function of *ξ* can be added to *M* without having any effect. Making use of Equation (95), Equation (104) simplifies to:

$$2\frac{\partial}{\partial\overline{\xi}^{\overline{\eta}}}\left[M+2i\nu\frac{\partial\Psi}{\partial\overline{\xi}^{\overline{\eta}}}\right]+\frac{u^{2}}{2}=0\,.\tag{105}$$

Finally, a second complex potential *χ* is introduced via:

$$M + 2\text{i}\nu \frac{\partial \Psi}{\partial \overline{\xi}^{\overline{\chi}}} = 2 \frac{\partial \chi}{\partial \overline{\xi}^{\overline{\chi}}} \, , \tag{106}$$

by which the two complex equations together with Equations (102) and (105) take the final form:

$$\frac{\partial \bar{u}\mu}{2} + \frac{p}{\varrho} + \mathcal{U} - \mathrm{i}\left[\frac{\partial \Psi}{\partial t} - 4\nu \frac{\partial^2 \Psi}{\partial \overline{\xi}^2 \partial \xi^t}\right] = 4 \frac{\partial^2 \chi}{\partial \overline{\xi}^2 \partial \xi^t} \tag{107}$$

$$-\frac{\mu^2}{2} = 4\frac{\partial^2 \chi}{\partial \overline{\xi}^2}.\tag{108}$$

By Equations (107) and (108), two complex equations are given containing the complex potential *χ*, the real-valued stream function Ψ and the pressure *p* as unknowns. They constitute a *first integral of Navier–Stokes equations* [73] in the sense that by taking the difference of the derivative of Equation (108) with respect to *ξ* with the derivative of Equation (107) with respect to *ξ*, the complex NS Equation (99) is recovered.

Compared to the classical complex variable approach outlined in Section 3.1, a complex valued potential field *χ* is required for the integration of the full 2D-NS equation in place of the real valued Airy stress function Φ utilised for the Stokes equation. A relationship between *χ* and Φ can be established for the particular case of *steady flow*: by setting *∂*Ψ/*∂t* = 0, Equation (107) simplifies to:

$$\frac{\bar{u}u}{2} + \frac{p}{\varrho} + \mathcal{U} = 4\frac{\partial^2}{\partial \overline{\varrho}^2 \partial \overline{\varrho}} \left[ \chi - \text{i}\nu \Psi \right] \text{ .}$$

which, apart from the additional nonlinear term *uu*¯ /2, corresponds to Equation (92) if Φ = *χ* − i*ν*Ψ is identified as the Airy stress function. Thus, for the steady flow case the two field Equations (107) and (108) simplify to:

$$\frac{\partial \bar{u}}{2} + \frac{p}{\varrho} + \mathcal{U} = 4 \frac{\partial^2 \Phi}{\partial \bar{\xi}^2 \partial \bar{\xi}'} \, \, \, \, \tag{109}$$

$$4\text{i}\nu \frac{\partial^2 \Psi}{\partial \overline{\xi}^2} + \frac{u^2}{2} = -4 \frac{\partial^2 \Phi}{\partial \overline{\xi}^2} \,'\,\tag{110}$$

in accordance with [38]. Note that Equation (109) takes the form of Bernoulli's equation apart from the second order derivative of the potential on its right hand side.

#### *3.3. Integration of the Dynamic Boundary Condition*

The mathematical derivation is completed by the specification of appropriate boundary conditions, which take the form of no-slip/no-penetration conditions at solid walls, inflow and outflow conditions and, in the case of film or multiphase flows, kinematic and dynamic boundary conditions at a free surface or internal interface. These are discussed in detail in [38,73]. However, a key feature of the above approach that requires emphasising is that the dynamic boundary condition associated with a free surface or internal interface can be similarly integrated in the case of steady flow leading to a corresponding first integral form that greatly simplifies enforcing such a condition compared to the standard method of addressing such problems employing the standard NS equations and

*Water* **2020**, *12*, 1241

boundary conditions written in primitive (i.e., observable) variables. Accordingly, this case as outlined briefly below.

Consider a simply connected domain with a boundary *ξ* = *f* (*s*), parametrised with respect to the arc length *s* of the boundary; furthermore, tangential and normal unit vectors along the boundary are taken to be *f* (*s*) and *n* = i*f* (*s*), respectively. In terms of the stream function Ψ, the kinematic boundary condition (*<sup>u</sup>* ¯ *f* ) = 0 ensures that along a free surface/interface the velocity in the normal direction vanishes, resulting in:

$$\frac{\partial \Psi}{\partial s} = 0.\tag{111}$$

The original dynamic boundary condition, after transformation into complex representation, reads:

$$\frac{p\_0 - p}{\varrho} n - 2i\nu \frac{\partial u}{\partial \overline{\xi}} \overline{f}' = \frac{\sigma}{\varrho} f''. \tag{112}$$

By making use of Equations (109) and (110), the pressure *p* and the viscous stresses i*ν∂u*/*∂ξ* can be replaced by derivatives of Φ, implying:

$$-\mu \frac{d\Psi}{\mathrm{d}s} + \left(\mathrm{d}l + \frac{p\_0}{\varrho}\right)n = \frac{\mathrm{d}}{\mathrm{d}s} \left[\frac{\sigma}{\varrho}f' + 4\mathrm{i}\frac{\partial\Phi}{\partial\overline{\_{\!\!\!}^{\!\!T}}}\right],\tag{113}$$

as the boundary condition for Φ. Since the specific potential energy *U* can be gauged with a constant, *U* + *p*0/ can be replaced by *U*. Making use of the kinematic boundary condition in Equation (111), integration of Equation (113) leads to:

$$4i\frac{\partial\Phi}{\partial\overline{\xi}} = -\frac{\sigma}{\varrho}f' + \int \mathcal{U}n\text{ds}\,\,\iota\tag{114}$$

as a first integral of the dynamic boundary condition. In [38], the same boundary condition is derived in real-valued form and where it is demonstrated that it can be decomposed into a Dirichlet and a Neumann boundary condition for Φ. The reduction of the original non-standard boundary condition in Equation (112) containing mixed contributions from different fields *u* and *p* to a mathematical standard form has proven to be extremely beneficial for the implementation of numerical methods of solution, as demonstrated in [6,38,74].

#### *3.4. Particular Flow Geometries as Exemplars*

Three different flow problems are used to demonstrate applications of the tensor approach described above; two of which lead to closed form analytical solutions, the third requiring a numerical solution.

#### 3.4.1. Uniaxial Flow: Flow over an Oscillating Plate

Consider the case of a uniaxial flow geometry with, *uy* = 0, or equivalently in complex formulation: *u* − *u*¯ = 0. With reference to Equation (95), the following PDE:

$$\frac{\partial \Psi}{\partial \xi} + \frac{\partial \Psi}{\partial \overline{\xi}} = 0$$

has to be satisfied, implying the following explicit forms:

$$\Psi^{\prime} = \,^{\prime}\Psi \left( \frac{\stackrel{\text{g}}{\xi} + \stackrel{\text{g}}{\xi}}{2\text{i}}, t \right) = \Psi^{\prime}(\mathcal{y}, t) \tag{115}$$

$$\mu \quad = \quad -2\mathbf{i}\frac{\partial \Psi}{\partial \overline{\xi}} = \Psi'(y, t) \tag{116}$$

for the stream function and the velocity; the prime denotes differentiation with respect to *y*. The general solution of Equation (108) is given by:

$$\chi = \mathbb{\mathcal{T}}w\_0(\xi, t) + w\_1(\xi, t) + \frac{1}{2} \int \left[ \int \Psi'(y, t)^2 \, \mathrm{d}y \right] \, \mathrm{d}y,\tag{117}$$

containing the two analytic functions *w*0(*ξ*, *t*), *w*1(*ξ*, *t*), being generalisations of the so-called *Goursat functions* [38]. Following the insertion of Equation (117), Equation (107) simplifies to:

$$\frac{dp}{d\xi} + lI - \text{i}\left[\dot{\Psi} - \nu \Psi''\right] = 4w\_0'(\xi, t) \text{ :} \tag{118}$$

the dot over a symbol denoting its time-derivative.

A horizontal plate of infinite extent covered by a fluid invokes a flow by a forced oscillatory movement, see Figure 2.

*U* cos(*ωt*)

**Figure 2.** Schematic showing the geometry for the flow over an oscillating plate.

The no-penetration condition is already fulfilled due to the prescribed flow configuration *uy* = 0, while the no-slip condition compels the fluid to mimic the oscillatory behaviour of the plate according to *ux*(*y*, 0) = *U* cos(*ωt*). By re-writing the latter in terms of the stream function Ψ via *ux* = *∂*Ψ/*∂y*, it takes the form of a Neumann boundary condition Ψ (0, *t*) = *U* cos(*ωt*). Since the stream function is gaugeable by an arbitrary constant, the additional Dirichlet condition Ψ(0, *t*) = 0 can be formulated without loss of generality. Additionally, the asymptotic condition Ψ (*y*, *t*) → 0 has to be fulfilled to ensure that the fluid tends to a state at rest far away from the plate.

Assuming vanishing pressure *p* = 0 and a wave-like solution of the form

$$\Psi = \mathcal{U}\Im \left\{ \exp(-\mathrm{i}\omega t) \frac{\exp(\mathrm{i}ky) - 1}{k} \right\} \,' \,. \tag{119}$$

together with fulfilment of the no-slip and no-penetration condition, then Equation (118) leads to the identity:

$$-\mathrm{i}lI\mathbb{S}\left\{\frac{-\mathrm{i}\omega+\nu k^{2}}{k}\exp(\mathrm{i}[ky-\omega t])+\mathrm{i}\frac{\omega}{k}\exp(-\mathrm{i}\omega t)\right\}=4w\_{0}'(\xi,t)',$$

requiring *k* and *ω* to satisfy the dispersion relation:

$$
\dot{\omega} = \nu k^2 \tag{120}
$$

for damped transverse waves in accordance to the classical result [22]. The Goursat function *w*<sup>0</sup> takes the form:

$$w\_0 = -\frac{\mathrm{i}}{4} \xi \nu l I \Im \left[ k \exp(-\mathrm{i}\omega t) \right]. \tag{121}$$

#### 3.4.2. Axisymmetric Flow: The Lamb-Oseen Vortex

Consider a class of flows given by the following particular form:

$$\Psi\_{\perp} = \Psi(\mathbf{r}, t) \,, \tag{122}$$

$$\mu \quad = \quad -2\mathbf{i}\frac{\partial \Psi}{\partial \overline{\pm}} = -\frac{\mathbf{i}\_{\sigma}^{\chi}}{r} \Psi'(r, t) \; , \tag{123}$$

$$
\sigma \quad := \quad \sqrt{\xi\_{\xi}^{\chi}},
\tag{124}
$$

of the stream function with the prime denoting differentiation with respect to *r*; that is flows, the streamlines of which form concentric circles, see Figure 3a. Assuming a particular solution for Equation (108) of the form *χ*<sup>p</sup> = *χ*p(*r*, *t*) leads to the the following expression:

$$r\frac{\partial}{\partial r}\left(\frac{\chi\_{\text{P}}'}{r}\right) = \frac{1}{2}\Psi'\left(r,t\right)^2,\tag{125}$$

and therefore:

$$\frac{\chi\_{\rm P}^{\prime}}{r} = \frac{1}{2} \int \frac{\Psi^{\prime}\left(r, t\right)^{2}}{r} d\mathbf{r} \,. \tag{126}$$

By inserting the above into Equation (107), the following simplified PDE results:

$$\frac{dp}{d\xi} + \mathcal{U} - \mathrm{i} \left[ \Psi - \frac{\nu}{r} \frac{\partial}{\partial r} \left( r \Psi' \right) \right] - \int \frac{\Psi' \left( r, t \right)^2}{r} \mathrm{d}r = 0 \,. \tag{127}$$

**Figure 3.** Vortex geometry (**a**) and time evolution of the velocity profile (**b**).

Via the introduction of the following similarity variable:

$$z = \frac{r}{\sqrt{\nu t}}$$

solutions of the form Ψ = *f*(*z*) are now searched for. Under the above assumptions the imaginary part of Equation (127) takes the form:

$$\frac{1}{t} \left[ \frac{z}{2} f'(z) + \frac{1}{z} \frac{\mathbf{d}}{\mathbf{d}z} \left( z f'(z) \right) \right] = 0$$

which, after making the substitution *g*(*z*) = *z f* (*z*), can be written conveniently as;

$$
\lg'(z) + \frac{z}{2}\lg(z) = 0,
$$

which, in turn, is easily solved by taking:

*<sup>g</sup>*(*z*) = *<sup>g</sup>*<sup>0</sup> exp −*z*2 4 , (128)

leading to:

$$f(z) = \text{g}\_0 \int \frac{1}{z} \exp\left(-\frac{z^2}{4}\right) \text{d}z = \frac{\text{g}\_0}{2} \text{Ei}\left(-\frac{z^2}{4}\right) \text{ ,} \tag{129}$$

where Ei denotes the integral exponential function. This particular solution contains a singularity; however, by considering the superposition:

$$\Psi = \frac{\mathcal{G}0}{2} \text{Ei}\left(-\frac{z^2}{4}\right) - \frac{\Gamma}{2\pi} \ln r\_{\nu}$$

with the well-known potential vortex <sup>Γ</sup> <sup>2</sup>*<sup>π</sup>* ln *r* as an alternative solution to Equation (127), the singularity is removed as follows: the complex velocity field resulting from Equation (123) reads

$$\mu = -\frac{\mathrm{i}\xi}{r}\Psi'(r,t) = \frac{\mathrm{i}}{r}\left[\frac{\Gamma}{2\pi} + \mathcal{g}\_0 \exp\left(-\frac{r^2}{4\nu t}\right)\right] \exp\left(\mathrm{i}\varphi\right) \omega$$

which is convergent in the limit *r* → 0 if and only if 2*πg*<sup>0</sup> = Γ. Finally, the solution reads:

$$u = \frac{\mathrm{i}\Gamma}{2\pi r} \left[1 - \exp\left(-\frac{r^2}{4\nu t}\right)\right] \exp\left(\mathrm{i}\rho\right) \,,\tag{130}$$

which is a reproduction of the classical Lamb-Oseen vortex [19]. The velocity profile |*u*| is shown in Figure 3b.

#### 3.4.3. Steady Film Flow over Topography

In the context of the numerical solution of flow problems, possibly involving a free surface, for which inertial effects cannot be ignored, the LSFEM, Cassidy [75], Thatcher [76], Bolton and Thatcher [77], has gained wide acceptability as a very effective and flexible approach—in particular when simple equal order elements in conjunction with highly efficient multigrid solvers are employed [6], exploiting the symmetry and positive definiteness of the resulting system matrices [78].

For reasons of practicality, the complex field Equations (109) and (110) are expressed as real-valued ones in terms of the velocity variables *ux* = *u* = *∂y*Ψ, *uy* = *u* = −*∂x*Ψ and first order derivatives, *φ<sup>x</sup>* = *∂x*Φ, *φ<sup>y</sup>* = *∂y*Φ, of the Airy stress function leading, together with the continuity Equation (53) and the condition:

$$\frac{\partial \phi\_1}{\partial y} - \frac{\partial \phi\_2}{\partial x} = 0,\tag{131}$$

to a system of four equations involving first order derivatives of *ux*, *uy*, *φx*, *φ<sup>y</sup>* only, conforming ideally to a first-order system least squares methodology.

Figure 4 shows the problem of gravity-driven film flow down a corrugated rigid surface inclined at an angle *α* to the horizontal considered by Marner [6]. Along the stationary corrugated surface, velocity Dirichlet conditions are imposed. While at the free surface, in addition to a kinematic boundary condition, two dynamic conditions are imposed resulting as inhomogeneous Dirichlet conditions for *φ*<sup>1</sup> and *φ*<sup>2</sup> from Equation (114) by decomposition into real and imaginary parts; these depend on the surface tension, the curvature and the potential energy density. The resulting free surface profile is

obtained by iterating over the kinematic condition while solving a sequence of flow problems with prescribed dynamic conditions in a fixed domain using the LSFEM.

**Figure 4.** Schematic of gravity-driven film flow down a corrugated rigid surface inclined at an angle *α* to the horizontal. The indices *i* and *j* in the boundary conditions run from 1 to 2.

Two representative results for two differently contoured and repeating surface shapes are shown in Figure 5:

**Figure 5.** Streamlines elucidating steady film flow over two different periodically corrugated inclined surfaces: on the left a smoothly varying feature, on the right a feature involving sharp, step changes.

The resulting streamline patterns reproduce exactly the observed corresponding experimental results of [79], obtained using silicon oil *Elbesil 145*: = 964.8 kg/m3, *ν* = 144.2 mm2/s, *σ* = 20.1 mN/m, *λ* = 20 mm, *A* = 4 mm, *H*<sup>0</sup> = 5 mm, *α* = 10◦.

Overall, the numerical method was shown [6] to produce an accurate and reliable result over a wide range of Reynolds and capillary numbers for the above problem.

#### *3.5. Tensor Potential Approach*

The obvious limitation of the 2D complex-variable formulation is its extension beyond two dimensions, since for the case of 3D viscous flow, a corresponding complex first integral formulation is, by definition, not possible; however, a real-valued tensor form is. In two dimensions, such a form can be established by decomposing both Equations (108) and (107) into real and imaginary parts and taking the linear combinations (107) ± (108), leading to four real-valued equations:

$$
\mu\_1^2 + \frac{p}{\varrho} + \mathcal{U} = 2\partial\_2^2 \Re \chi + 2\partial\_1 \partial\_2 \Im \chi \tag{132}
$$

$$
\mu\_2^2 + \frac{p}{\varrho} + \mathcal{U} = 2\partial\_1^2 \Re \chi - 2\partial\_1 \partial\_2 \Im \chi \,\, \, \, \, \tag{133}
$$

$$
\mu\_1 \mu\_2 = -\left\{ \partial\_1^2 - \partial\_2^2 \right\} \odot \chi - 2\partial\_1 \partial\_2 \Re \chi \,\tag{134}
$$

$$
\partial\_t \Psi = -\nabla^2 \left[ \Im \chi - \nu \Psi \right],
\tag{135}
$$

*Water* **2020**, *12*, 1241

which, in the context of extension to higher dimensions, can be written in the following convenient and compact tensor form:

$$
\mu\_i \mu\_j + \left[\frac{p}{\varrho} + \mathcal{U}\right] \delta\_{ij} = -\partial\_k \partial\_k \overline{a}\_{ij} + \partial\_i A\_j + \partial\_j A\_i - \partial\_k A\_k \delta\_{ij},\tag{136}
$$

with the tensor and vector fields given by *a*˜*ij* = −*χδij* and *Aj* = *∂ia*˜*ij* + *εjk∂kχ*, respectively. The abbreviations *∂<sup>k</sup>* := *∂*/*∂xk* and *∂<sup>t</sup>* := *∂*/*∂t* imply partial differentiation, while *δij* and *εij* denote the Kronecker delta function and the two-dimensional Levi-Civita symbol, respectively.

Recently [39], the corresponding tensor potential formulation for three dimensions has been established, an essential underpinning being analogies drawn with the methodical reduction of Maxwell's equations [6]. The continuity Equation (53) is shown to be fulfilled identically following the introduction a streamfunction vector Ψ*<sup>k</sup>* for the velocity, in accordance with *ui* = *εijk∂j*Ψ*<sup>k</sup>* and *εijk* denoting the corresponding 3D Levi-Civita symbol, representing a 3D generalisation of the 2D streamfunction [21]. Compared to the 2D first integral of the NS equations, the corresponding 3D formulation utilises an independent symmetric tensor potential *a*˜*ij* and a vector potential *ϕ<sup>i</sup>* with the indices *i*, *j* taking values from 1 to 3. The auxiliary vector field *Aj* in the 3D case reads:

$$A\_j = \partial\_k \vec{a}\_{kj} + \varepsilon\_{jlk} \partial\_l \varphi\_k \,. \tag{137}$$

and Equation (135) has to replaced by its corresponding 3D form:

$$
\partial\_l \Psi\_{\,\,l} = \nu \partial\_k \partial\_k \Psi\_{\,\,l} - \varepsilon\_{nkl} \partial\_k \partial\_{\,m} \overline{a}\_{ml} \,. \tag{138}
$$

Like the Clebsch variables considered in Section 2, the tensor potential *a*˜*ij* and the vector potential *ϕ<sup>i</sup>* are not unique and can be gauged in a beneficial way. Obviously, by performing the operations:

$$\vec{a}\_{i\dot{j}} \longrightarrow \vec{a}\_{i\dot{j}} + \vec{\partial}\_{\dot{i}} \vec{a}\_{\dot{j}} + \vec{\partial}\_{\dot{j}} \vec{a}\_{i\dot{i}} - \vec{\partial}\_{k} \vec{a}\_{k} \delta\_{\dot{i}\dot{j}},\tag{139}$$

$$
\varphi\_i \longrightarrow \varphi\_i + \partial\_i \zeta \,, \tag{140}
$$

for an arbitrary vector field *α<sup>i</sup>* and an arbitrary scalar field *ζ*, the field Equations (136) remain invariant. These rules are utilised in [39] to establish bona fide gauging scenarios. In the same paper the prescription of appropriate commonly occurring physical and necessary auxiliary boundary conditions, incorporating for completeness the derivation of a first integral of the dynamic boundary condition at a free surface, is established, together with how the general approach can be advantageously reformulated for application in solving unsteady flow problems with periodic boundaries.

Using a tensor formulation, the approach is suitable for use in the case of an arbitrary number of dimensions: in [80], a potential-based first integral form is established for the 4D energy-momentum equations for flows under relativistic conditions.

#### **4. Discussion**

Based on a detailed analysis and discourse, the two different approaches considered above can be explained in the light of their different origins: the Clebsch representation of the velocity can, according to recent analysis [56], be understood as a natural outcome of Galilean invariance; whereas, Airy's stress function originates historically from the 2D static equilibrium of internal forces. In the course of a long and growing series of research contributions both approaches have been generalised and made available for use in solving arbitrary flow problems. A Clebsch transformation has emerged that applies to arbitrary forces in the equations of motion, including viscous ones, while extensions to the Airy stress function approach applies to cases beyond static equilibrium, in particular to unsteady flows with inertia and, after rearrangement of the complex equations, to tensor equations (in terms of a tensor potential) for 3D viscous flows and even for the case of 4D relativistic flows. The use of a tensor potential has parallels to Maxwell's theory of electromagnetism [6,39,80]. The capabilities of both approaches have been convincingly demonstrated through the solution of a variety of illustrative examples.

Despite the very positive stage of development of both methods, some open questions remain: the first is that the field Equations (64)–(68) resulting, via the extended Clebsch transformation according to Theorem 2, from the NS equations are not self-adjoint; as in the inviscid flow case. The non-existence of a Lagrangian appears to be due to energy dissipation caused by viscosity. Anthony [81] poses a possible strategy to overcome this problem by including thermal degrees of freedom and related inner energy in order to remain consistent with Noether's theorem, which implies conservation of energy for Lagrangians being invariant with respect to time-translations. Note that the present work, by *dissipation*, the irreversible transfer of energy from mechanical to thermal energy is understood from the physicist's viewpoint, while the total energy (the sum of mechanical and thermal energy) is conserved. Based on Seliger and Witham [24]'s classical work, Zuckerwar and Ash [82,83] suggest a Lagrangian considering only volume viscosity, leading to equations of motion containing, only qualitatively, the effect of volume viscosity but differing quantitatively from the compressible NS equations—also known as the Navier–Stokes–Duhem equations [84,85] without shear viscosity. They interpret their result as a generalisation of the theory of viscous flow towards thermodynamic non-equilibrium. A Lagrangian for viscous flow considering both shear viscosity and volume viscosity has been suggested by Scholle and Marner [86], Marner et al. [87], again leading to equations of motion that differ from the Navier–Stokes–Duhem equations by non-equilibrium terms. A striking feature of their Lagrangian is a *discontinuity*, causing fluctuations on a microscopic scale and revealing parallels to a stochastic variational description as in statistical physics; see, for example, [88–93]. Since these considerations go beyond the scope of classical fluid mechanics, further research on this particular field is required.

A second unanswered question is whether a general and all-encompassing potential approach exists reducing to both the Clebsch and the tensor potential approach as special cases. The search for this 'missing link' between two conceptually different approaches represents another future research topic of general interest. A promising first step could be an analysis of the Clebsch transformation in the sense of general relativity [94], followed by a particular classical limit as shown in Lightman et al. [95], problem 5.31 pp. 35, 227–228.

Thirdly, the first integral of the 4D energy-momentum equations based on a tensor potential [80] points to the future use of mathematical techniques and methods of solutions not currently applicable to the field equations in their original form, in particular the use of matrix structures within the framework of Clifford algebra, based on quaternions or Dirac matrices with the goal of developing highly efficient methods of solution. Having mapped the entire problem to a matrix-algebra framework, the limit *c* → ∞ could be applied in order to provide efficient solutions of the classical NS equations. Although implementation of such a matrix-algebra techniques remains speculative at this stage, it deserves further investigation since its utilisation can lead to significant economic gains in the computation of fluid flows, in a similar fashion to the use of quaternions representing spatial rotation operations [96,97], potentially leading to the formulation of highly efficient and predictive CFD software.

**Author Contributions:** Conceptualisation, M.S. and P.H.G.; methodology, M.S. and F.M.; software, F.M.; investigation, F.M. and P.H.G.; writing—original draft preparation, M.S.; writing—review and editing, P.H.G. and F.M.; visualisation, M.S. and F.M.; funding acquisition, M.S., P.H.G. and F.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the German Research Foundation (DFG) grant numbers SCHO 767/6-1 and SCHO 767/6-3.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A. Proof of the Existence of a Representation with Two Pairs of Clebsch Variables**

Lin [44,46] proposed a variational principle for superfluid helium resulting by variation in a Clebsch representation with three pairs of variables:

$$\vec{u} = \nabla \Phi + \lambda\_1 \nabla \xi\_1 + \lambda\_2 \nabla \xi\_2 + \lambda\_3 \nabla \xi\_3,\tag{A1}$$

where the three fields *ξ<sup>i</sup>* (*i* = 1, 2, 3) are identified physically as material coordinates while the three conjugated fields *λ<sup>i</sup>* have been introduced as Lagrange multipliers. It can be proven easily that for a fluid flow a representation of the velocity in the form of Equation (A1) always exist, even if Φ = 0 is assumed. The essence is the existence of both a field representation *ψ* = *ψ* (*xi*, *t*) and material representation *ψ* = Ψ *ξj*, *t* for an arbitrary field *ψ* and a related invertible transformation:

$$\mathbf{x}\_{i} = \mathbf{x}\_{i} \left( \mathbf{f}\_{j}, \mathbf{t} \right), \tag{A2}$$

between both [22,98]. By taking the gradient of Equation (A2), the relation:

$$
\vec{\varepsilon}\_{i} = \nabla x\_{i} = \frac{\partial x\_{i}}{\partial \vec{\xi}\_{j}} \nabla \vec{\xi}\_{j} \tag{A3}
$$

is obtained with the deformation gradient *<sup>∂</sup>xi ∂ξ<sup>j</sup>* . Thus for an arbitrary field *u* = *ui ei* it follows that:

$$\vec{u} = u\_i \vec{e}\_i = \underbrace{u\_i \frac{\partial x\_i}{\partial \xi\_j^x}}\_{\lambda\_j} \nabla \xi\_{j\_{\prime}} \tag{A4}$$

corresponding to the form (A1) with Φ = 0.

Next the material representation *λ*<sup>3</sup> = Λ<sup>3</sup> *ξj*, *t* for the variable *λ*<sup>3</sup> is utilised in order to derive the identity:

$$\begin{split} \vec{u} &= \lambda\_1 \nabla\_{\xi}^{\mathfrak{x}} + \lambda\_2 \nabla\_{\xi}^{\mathfrak{x}} + \Lambda\_3 \left( \mathfrak{z}\_j^{\mathfrak{x}}, t \right) \nabla\_{\xi}^{\mathfrak{x}} \\ &= \nabla \int \Lambda\_3 \left( \mathfrak{z}\_j^{\mathfrak{x}}, t \right) \mathrm{d}\_{\xi}^{\mathfrak{x}} + \left[ \lambda\_1 - \frac{\partial}{\partial \xi\_1} \int \Lambda\_3 \left( \mathfrak{z}\_j^{\mathfrak{x}}, t \right) \mathrm{d}\_{\xi}^{\mathfrak{x}} \right] \nabla \xi\_1^{\mathfrak{x}} + \left[ \lambda\_2 - \frac{\partial}{\partial \xi\_2} \int \Lambda\_3 \left( \mathfrak{z}\_j, t \right) \mathrm{d}\_{\xi}^{\mathfrak{x}} \right] \nabla \xi\_2^{\mathfrak{x}}, \end{split}$$

corresponding via the definitions:

$$\mathfrak{gl} := \int \Lambda\_3 \left( \mathfrak{z}\_{j\_\prime} t \right) \operatorname{d} \mathfrak{z}\_{\mathfrak{z}} \Big|\_{\mathfrak{z}\_j = \mathfrak{z}\_j(x\_i, t)} \tag{A5}$$

$$\mathfrak{a}\_1 := \lambda\_1 - \frac{\mathfrak{d}}{\partial \mathfrak{J}\_1} \int \Lambda\_3 \left( \mathfrak{J}\_j^{\mathfrak{x}}, t \right) \operatorname{d} \mathfrak{J}\_3 \Big|\_{\mathfrak{J}\_j^{\mathfrak{x}} = \mathfrak{J}\_j^{\mathfrak{x}}(x\_i, t)} \tag{A6}$$

$$\mathcal{A}\_2 := \lambda\_2 - \frac{\partial}{\partial \xi\_2^x} \int \Lambda\_3 \left( \xi\_j^x, t \right) \mathrm{d}\xi\_3^x \Big|\_{\xi\_j = \xi\_j^x(x\_i, t)} \tag{A7}$$

$$\mathcal{B}\_1 := \mathcal{J}\_1 \tag{A8}$$

$$\beta\_2 := \S\_2 \tag{A9}$$

*Water* **2020**, *12*, 1241

to a Clebsch representation with two pairs of Clebsch variables:

$$\vec{u} = \nabla \varphi + \mathfrak{a}\_1 \nabla \beta\_1 + \mathfrak{a}\_2 \nabla \beta\_2 \tag{A10}$$

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Stochastic Approaches to Deterministic Fluid Dynamics: A Selective Review**

#### **Ana Bela Cruzeiro**

Departamento Matemática I.S.T. and Grupo de Física-Matemática University, Av. Rovisco Pais, 1049-001 Lisbon, Portugal; ana.cruzeiro@tecnico.ulisboa.pt

Received: 20 January 2020; Accepted: 13 March 2020; Published: 19 March 2020

**Abstract:** We present a stochastic Lagrangian view of fluid dynamics. The velocity solving the deterministic Navier–Stokes equation is regarded as a mean time derivative taken over stochastic Lagrangian paths and the equations of motion are critical points of an associated stochastic action functional involving the kinetic energy computed over random paths. Thus the deterministic Navier–Stokes equation is obtained via a variational principle. The pressure can be regarded as a Lagrange multiplier. The approach is based on Itô's stochastic calculus. Different related probabilistic methods to study the Navier–Stokes equation are discussed. We also consider Navier–Stokes equations perturbed by random terms, which we derive by means of a variational principle.

**Keywords:** Navier–Stokes equation; stochastic Lagrangian flows; stochastic variational principles; stochastic geometric mechanics

#### **1. Introduction**

The dynamics of an incompressible viscous fluid is modeled by the Navier–Stokes equation, a second order, nonlinear, partial differential equation describing the balance of mass and momentum of the fluid flow. In the absence of external forces and considering a perfectly incompressible fluid, the Navier–Stokes equation reads

$$\frac{\partial}{\partial t}u + (u, \nabla)u = \nu \Delta u - \nabla p \tag{1}$$

where the velocity field u is required to satisfy the incompressibility condition div*u* = 0, the fluid density being equal to 1. The constant *ν* denotes the kinematic viscosity and *u* the fluid velocity. Moreover, the symbol *p* stands for the pressure within the fluid, and is yet another unknown in the equation. Equation (1) has a huge number of applications in physics and engineering.

The Navier–Stokes equation is deterministic, but it is well known that some of its solutions seem to exhibit random behavior, which might eventually provide an insight into the onset of turbulence although from a purely mathematical point of view, the very important problem of existence and smoothness of the solutions to Equation (1) remains largely unsolved to this day. In this article we will review various ways of introducing randomness into the analysis of the solutions to Equation (1).

In Physics, the most famous deterministic equation hiding randomness is the Schrödinger equation. There is no mathematical probability theory behind it, but there is manifestly randomness in any Quantum Physics Lab. This puzzling situation did not prevent Quantum Theory to develop an impressive corpus of techniques to control, in particular, the transition from (in principle differentiable) solutions of classical equations of motion to some forms of randomness. Is it possible to draw an analogy between this classical/quantum relation and the one between Euler and Navier–Stokes equations? Of course the status of the two hydrodynamical equations are quite different from the

above ones. Euler equation, although modeling only "dry water" (Von Neumann) is already quite complicated. It does not seem, yet, to have proof that it does not produce singularities. But the abovementioned analogy could, in particular, provide an insight into the onset of turbulence.

There are different ways to introduce randomness: uncertainty may have its origin in errors in the initial conditions, for example. In this case statistical approaches are considered: one studies the time evolution of some probability measure, supported on the relevant physical initial data. This is part of the statistical approach to turbulence, initiated already in the 19th century (c.f, among many others, [1,2]). Other various types of Langevin dynamics, using stochastic diffusion processes, have been proposed to describe equilibrium and non-equilibrium dynamics as well as Kraichnan's model in turbulent advection (c.f, for instance, [3]). On the other hand uncertainty may be generated by the chosen numerical model (we refer to [4] for a discussion on these issues in climate modeling, see also [5]).

Another very popular way to introduce stochasticity is to perturb the Navier–Stokes equation with random forces, so that stochastic partial differential equations then come into play. There is a huge literature in this subject, after the pioneering mathematical work [6]. Stochasticity is typically introduced at a Eulerian level, although some stochastic Lagrangian models of Langevin type (with smooth Lagrangian trajectories and stochastic velocities) have also been considered in turbulence ([7]).

More recently stochastic advection by Lie transport was introduced by D.D. Holm in [8]. The resulting equations of motion are stochastic partial differential equations and the approach is also Eulerian.

It would be impossible to mention here all stochastic approaches to fluid dynamics. We have only chosen some topics and a few corresponding references. It is worthwhile to mention that there are some interesting probabilistic representation formulae. Representing solutions of partial differential equations as expected values of functionals of stochastic processes is a tradition in the field of stochastic analysis and has also been considered for fluid dynamics. Very roughly speaking, we can find three different approaches: the probabilistic representation of the vorticity field as in [9], the analysis through branching processes and the Fourier transform as in [10] and that using Lagrangian diffusion processes as in [11].

On the other hand, in this paper, we are concerned with variational principles for deterministic dissipative fluid dynamics, and a brief state-of-the-art review follows. There are relatively few references on the special stochastic view of deterministic fluid dynamics advocated here. The oldest we know are [12,13], where the Laplacian term of the Navier–Stokes equation was interpreted, quite informally, as the presence of an underlying Brownian motion. In the famous paper [14], as well as in [15,16], a rigorous geometric strategy for the Euler equation was developed. These geometric ideas have been extended by us to the Navier–Stokes equation, together with the associated variational principle, giving rise to a new stochastic geometric approach to dissipative dynamics. Recently, [17,18], more physical-oriented works, were influenced by references [12,13] and also by our work.

The review paper presented here is, of course, very selective, not in terms of the qualities of quoted references, but in terms of their perspectives, in relation with the interplay between deterministic and stochastic viewpoints. We apologize to the many authors whose important works we were not able to describe here. We hope, however, to describe a special viewpoint in a consistent manner.

In a nutshell, the goal of this review is to show that probabilistic methods play a central role in the study of the deterministic Navier–Stokes equation, both from a conceptual and a practical point of view. Indeed, on the one hand the solutions to that equation satisfy stochastic variational principles. On the other hand, since the solutions to the Navier–Stokes equation may be interpreted as drifts of diffusion processes, one may import many techniques from stochastic analysis to investigate them. Those techniques include the use of stochastic differential equations, of forward–backward stochastic systems and of related numerical techniques.

#### **2. Results**

As a dissipative system, there are well known obstructions if we want to derive the classic Navier–Stokes equation from a deterministic variational principle. Allowing the Lagrangian paths to be random and using Itô differential calculus (c.f [19]), we show how we can still derive the Navier–Stokes equation, without any randomness in the external forces, from a (stochastic) variational principle, where the Lagrangian functional is still the classic one, but computed over random Lagrangian trajectories. In some sense, the equation can be regarded as a generalized geodesic equation defined in some suitable space and the variational principle reduces to Hamilton's principle for the incompressible Euler equation when the viscosity vanishes. Our variational approach can be extended in order to consider stochastic Navier–Stokes equations as well.

After recalling Arnold's variational approach to the Euler equation, we describe in Section 4 two stochastic variational principles for the Navier–Stokes equation. The first is a "direct" generalization of Arnold's variational principle to the viscous case, the viscosity being associated with the random behavior of the particles. Incompressibility is incorporated in the definition of their trajectories. The second imposes incompressibility via a Lagrange multiplier.

Having justified our approach of Navier–Stokes equation using stochastic Lagrangian paths, Section 5 is devoted to several possible mathematical methods to study such paths. We can use the theory of forward–backward stochastic differential equations: this is explained in Section 5.1. We can also use entropy methods, since our action functional is essentially given by an entropy quantity. A notion of weak solutions, in the spirit of Brenier's work for the Euler equation and of optimal transport theory, allows us to consider cases where other methods are not accessible by lack of regularity. In Section 6 we describe some stability properties of the Navier–Stokes stochastic Lagrangian flows. The following paragraph is devoted to stochastic perturbations of the Navier–Stokes equation: we show that they can also be derived from a variational principle. Other equations and methods are mentioned in Section 8, as well as some future research problems. Finally a brief appendix contains basic notions of Itô stochastic calculus that are used in this paper.

#### **3. The Non Viscous Case**

As pointed out by Arnold ([14]), the incompressible Euler equation, corresponding to the case where the viscosity of the fluid is equal to zero, is an equation of geodesics on a suitable (infinite dimensional) space of functions. If we consider Lagrangian paths *g*(*t*, *x*) such that *<sup>∂</sup> <sup>∂</sup><sup>t</sup> g*(*t*, *x*) = *u*(*t*, *g*(*t*, *x*)), with *u* solution of the Euler equation

$$\frac{\partial}{\partial t}u + (u.\nabla)u = -\nabla p\_\prime \quad \text{div}\, u = 0,\tag{2}$$

the acceleration *<sup>∂</sup>*<sup>2</sup> *<sup>∂</sup><sup>t</sup> <sup>g</sup>* is equal to a gradient function and thus, at every time, orthogonal (for the *<sup>L</sup>*<sup>2</sup> scalar product) to vector fields of zero divergence. This means that, if we endow the space of diffeomorphisms preserving the volume measure *m* of the underlying configuration space a structure of manifold, the Lagrangian flows *g*(*t*, ·) will be geodesics in such a manifold since its tangent space will consist of divergence free vector fields. In particular they will be critical points of the action functional defined by the kinetical energy:

$$S[\mathcal{g}] = \frac{1}{2} \int\_0^T \int |\dot{\mathcal{g}}(t, \mathbf{x})|^2 dm(\mathbf{x}) dt = \frac{1}{2} \int\_0^T \left\| |\dot{\mathcal{g}}(t)| \right\|\_{L^2(dm)}^2 dt. \tag{3}$$

More precisely, Arnold showed that the Euler equation above corresponds to the equation of the geodesic flow of the (right-invariant) *L*<sup>2</sup> metric on the group of diffeomorphisms of the underlying configuration space that preserve the volume measure and have a certain Sobolev regularity. The velocity can be recovered from the Lagrangian flow: *u*(*t*, *x*)=( *<sup>∂</sup> <sup>∂</sup><sup>t</sup> gt*)(*g*−<sup>1</sup> *<sup>t</sup>* (*x*)).

This program was rigorously developed in [15] and geodesics were shown to exist locally and under certain regularity restrictions on the initial conditions. Moreover the information on the geometry of the problem had important consequences, for instance in describing the chaotic behavior of Euler Equation (c.f [16]) and its consequences for weather prediction. Instability of the geodesic Euler flows can be described in terms of the sectional curvatures of the group of volume preserving diffeomorphisms. Explicit estimates for the curvatures, being non positive, show that, essentially, the weather is unpredictable.

Notice that this is a special case of Lagrangian system treated in Geometric Mechanics via variational principles on Lie groups ([20]). Indeed the space of volume preserving diffeomorphisms has also a group structure for the composition of maps and the *L*<sup>2</sup> metric is right-invariant.

#### **4. The Deterministic Navier–Stokes Equation**

The situation concerning variational principles for the Navier–Stokes equation is radically different since, contrary to the Euler one, this system is dissipative. Our approach describes dissipation through an introduction of noise in the Lagrangian trajectories. We replace deterministic Lagrangian paths by semimartingales of the form

$$d\mathfrak{J}\_l^x(\mathbf{x}) = dM\_l(\mathbf{x}) + D\_l \mathfrak{J}\_l(\mathbf{x}), \qquad \mathfrak{J}\_0(\mathbf{x}) = \mathbf{x} \tag{4}$$

where *Mt* is a martingale (for instance, a Brownian motion) and *Dtξ<sup>t</sup>* denotes the bounded variation part of the semimartingale. The use of the notation *Dt* is not an accident, since it corresponds to a mean derivative in time (recall that, almost-everywhere, *ξ<sup>t</sup>* satisfying Equation (3) will not be time-differentiable). Consider the definition of the generalized derivative *Dt*. Namely, for every regular function *F*,

$$D\_t F(t, \xi\_t^\tau) = \lim\_{\varepsilon \to 0} \frac{1}{\varepsilon} E\_t[F(t + \varepsilon, \xi\_{t+\varepsilon}^\tau) - F(t, \xi\_t^\tau)],\tag{5}$$

where *Et* denotes the conditional expectation with respect to the past information of the process *ξ*. Applying this definition to the identity function, and since the operator *Dt* vanishes when considered on martingales, we see that *Dtξ<sup>t</sup>* coincides with the drift of the process and can indeed be understood as a regularized derivative. This notion has been known in Stochastic Analysis since its beginnings, as it corresponds to the definition of the generator of the process *ξt*, but it became more relevant in dynamics with the works of Nelson [21].

When *ξ* takes values in a non-Euclidean space, notably a Riemannian manifold or a Lie group, one can extend the definition of generalized derivative using parallel transport. In this paper we will mainly consider the Euclidean setting. We denote the configuration space by O and, for now, we assume that O has no boundary (case with periodic boundary conditions, for example).

#### *4.1. Stochastic Geometric Mechanics Approach*

We consider the space *GV* of bijective maps *<sup>g</sup>* : O→O which belong to *<sup>L</sup>*<sup>2</sup> = *<sup>L</sup>*2(*dm*) and keep the volume measure invariant, namely such that *f*(*g*(*x*))*dm*(*x*) = *f dm*(*x*) for every (continuous) function *f* . This infinite dimensional space of maps has a group structure under composition and a Riemannian one with respect to the *L*<sup>2</sup> metric ([15]). The tangent space to the identity map consists of vector fields on O with zero divergence. Moreover a simple covariant derivative can be defined: for such two vector fields *X* and *Y*, ∇*XY* = Π(*∂XY*) where *∂* is the usual Lie derivative and Π the projection operator from *L*<sup>2</sup> to divergence free vectors associated to the Helmholtz (or Hodge, according to the source) decomposition.

The next result was proved in [22], then extended to Riemannian manifolds in [23] and formulated on general Lie groups in [24].

We consider a stochastic action functional where the Lagrangian is the kinetic energy computed on the generalized derivative of a semimartingale. The corresponding norm is the *L*<sup>2</sup> one. More precisely, for a *GV*-valued semimartingale *ξ* as in Equation (4) we define

$$S[\mathcal{J}] = \frac{1}{2} E \int\_0^T \int\_{\mathcal{O}} |D\_t \mathbf{j}\_t(\mathbf{x})|^2 dt dm(\mathbf{x}) \tag{6}$$

where *E* means expectation with respect to the underlying probability. A particular class of such semimartingales is the following. To a time dependent vector field with zero divergence *u*(*t*, ·), *<sup>t</sup>* ∈ [0, *<sup>T</sup>*], and belonging to *<sup>L</sup>*2, we associate the stochastic differential equation (c.f [19], for example, as a reference for Itô stochastic calculus, as well as the Appendix A in this article),

$$d\mathbf{g}\_t^{\mu}(\mathbf{x}) = \sqrt{2\nu} \, d\mathcal{W}\_t + \mu(t, \mathbf{g}\_t^{\mu}(\mathbf{x}))dt,\qquad \mathbf{g}\_0^{\mu}(\mathbf{x}) = \mathbf{x} \tag{7}$$

where *W* is a standard Brownian motion. This equation defines a stochastic flow of maps on O that belong to *GV* and satisfy *Dtg<sup>u</sup> <sup>t</sup>* = *u*(*t*, *gt*).

Considering a simple Brownian motion may be regarded as oversimplifying. Actually one should introduce martingales driven by vector fields modeling the correlations observed in the physical model ([8]). Typically Lagrangian paths are of the form

$$d\mathcal{g}\_t = \sum\_k H\_k(\mathcal{g}\_t) d\mathcal{W}\_t^k + \mu(t, \mathcal{g}\_t) dt$$

where *Hk* are correlation eigenvectors and *W<sup>k</sup>* independent Brownian motions. Our approach covers such cases ([22,24,25]) and we chose here a Brownian motion only to simplify the exposition.

We are interested in derivating *S*[*gu*]; in particular we want to consider variations of the paths *g<sup>u</sup>* for which the functional above is still well defined, i.e., they are still *GV*-valued semimartingales. Consider the exponential type functions

$$e\_t(\epsilon v)(x) = x + \epsilon \int\_0^t \dot{v}(s, e\_s(\epsilon v)(x)) ds$$

with  > 0 and where *v*(*t*, ·) is a smooth time dependent vector field such that *v*(0) = *v*(*T*) = 0 and div *v*(*t*, ·) = 0 for every *t* ∈ [0, *T*]. Notice that, up to the first order in , we have *et*( *v*)(*x*) *x* +  *v*(*t*, *x*). The variations of the paths *gu*(*t*) will be defined by left composition, since the functional is right-invariant:

$$\mathbf{g}\_t^{\mathfrak{u},\mathfrak{c}} = e\_\mathfrak{l}(\boldsymbol{\epsilon}\boldsymbol{\upsilon}) \circ \boldsymbol{\mathcal{g}}^\mathfrak{u}(t),$$

We have, using Itô calculus,

$$d\mathcal{g}\_t^{\boldsymbol{\mu},\boldsymbol{\varepsilon}} = \nabla e\_l(\boldsymbol{\varepsilon}\boldsymbol{\upsilon})(\boldsymbol{g}\_t^{\boldsymbol{\mu}})\sqrt{2\boldsymbol{\nu}}d\mathcal{W}\_l + [\dot{e}\_l(\boldsymbol{\varepsilon}\boldsymbol{\upsilon}) + (\boldsymbol{\mu}.\nabla)e\_l(\boldsymbol{\varepsilon}\boldsymbol{\upsilon}) + \boldsymbol{\nu}\Delta e\_l(\boldsymbol{\varepsilon}\boldsymbol{\upsilon})](\boldsymbol{g}\_t^{\boldsymbol{\mu}})]dt$$

By the definition of *et*(*v*),

$$\frac{d}{d\epsilon}\Big|\_{\epsilon=0} S[e, (\epsilon v) \circ \mathcal{g}^u(\cdot)] = E \int\_0^T \left( \int D\_t \mathcal{g}^u(t)(\mathbf{x}) \, D\_t v(\mathcal{g}^u(t)(\mathbf{x})) dm(\mathbf{x}) \right) dt$$

and by Itô's formula,

$$\begin{aligned} \int d\int \mathcal{D}\_l g^u(t)(\mathbf{x}) . v(g^u(t)(\mathbf{x})) dm(\mathbf{x}) &= \int d\mathcal{D}g^u(t)(\mathbf{x}) . v(g^u u(t)(\mathbf{x})) dm(\mathbf{x}) + \int \mathcal{D}g^u(t)(\mathbf{x}) . dv(g^u(t)(\mathbf{x})) dm(\mathbf{x}) \\\\ &+ \int d\mathcal{D}\_l g^u(t)(\mathbf{x}) . dv(g^u(t)(\mathbf{x})) dm(\mathbf{x}) \end{aligned}$$

The last (Itô's contraction) term is equal to

$$2\nu(\int (\nabla v \otimes \nabla u) (g^u(t)(\mathfrak{x})) dm(\mathfrak{x})) dt$$

As *v*(0) = *v*(*T*) = 0 this implies,

$$\frac{d}{d\varepsilon}\Big|\_{\varepsilon=0} \mathcal{S}[e, (\varepsilon v) \circ \mathcal{g}^{\mu}(\cdot)] = -E \int\_0^T (\int (D\_t D\_t \mathcal{g}^{\mu}(t)(\mathbf{x}) dm(\mathbf{x})) dt - 2\nu E \int\_0^T (\int (\nabla v \otimes \nabla u)(\mathcal{g}^{\mu}(t)(\mathbf{x})) dm(\mathbf{x})) dt$$

On the other hand

$$D\_t D\_t \mathcal{g}^\mu(t) = D\_t \mu(t, \mathcal{g}^\mu(t)) = (\frac{\partial}{\partial t} \mu + (\mu . \nabla) \mu + \nu \Delta \mu)(\mathcal{g}^\mu(t)).$$

Therefore, using the invariance of the measure with respect to the process *g<sup>u</sup>* and integration by parts, we obtain

$$\frac{d}{d\varepsilon}\bigg|\_{\varepsilon=0} S[\varepsilon(\varepsilon v) \diamond g^u(\cdot)] = -E \int\_0^T \int \left( ([\frac{\partial}{\partial t}u + (u.\nabla)u - \nu \Delta u].v)(t, g^u(t)(x))dm(x) \right) dt$$

$$= -\int\_0^T \left( \int ([\frac{\partial}{\partial t}u + (u.\nabla)u - \nu \Delta u].v)(t,x)dm(x) \right) dt$$

for every *v* with zero divergence, which means that *<sup>∂</sup> <sup>∂</sup><sup>t</sup> u* + (*u*.∇)*u* − *ν*Δ*u* is the gradient of some function.

We have therefore the following

**Theorem 1.** *Let <sup>u</sup>*(*t*, ·) *be a smooth time dependent divergence-free vector field defined on* [0, *<sup>T</sup>*] × O*. Let <sup>g</sup>u*(*t*) *be a stochastic Brownian flow with diffusion coefficient* <sup>√</sup>2*<sup>ν</sup> and drift <sup>u</sup> (as in Equation (7)). Then <sup>g</sup><sup>u</sup> is a critical point of the energy functional S if and only if there exists a function p such that u*(*t*) *verifies the incompressible Navier–Stokes Equation (1).*

**Remark 1.** *When* O *is a Riemannian manifold (say, without boundary), we can still define an action functional of the form in Equation (6), but some more concepts are needed.*

In that case the Itô differential of an O-valued semimartingale *Y* is defined by

$$dY\_t = P\left(Y\right)\_t d\left(\int\_0^\cdot P\left(Y\right)\_s^{-1} \circ dY\_s\right)\_t$$

where *P* (*Y*)*<sup>t</sup>* : *TY*<sup>0</sup> *M* → *TYt M* is the parallel transport associated with the Levi–Civita connection along *t* → *Yt*. Alternatively, in local coordinates,

$$d\mathcal{Y}\_t = \left(d\mathcal{Y}\_t^i + \frac{1}{2}\Gamma^i\_{jk}(\mathcal{Y}\_t)d\mathcal{Y}\_t^j \otimes d\mathcal{Y}\_t^k\right)\partial\_i\mathcal{Y}\_t^j$$

where Γ*<sup>i</sup> jk* are the Christoffel symbols of this connection.

If the semimartingale *Yt* has an absolutely continuous drift, we denote it by *DYt*: for every 1-form *<sup>α</sup>* <sup>∈</sup> <sup>Γ</sup>(*T*∗O), the finite variation part of · <sup>0</sup> *α*(*Yt*), *dYt* is · <sup>0</sup> *α*(*Yt*), *DYt dt*. We consider an incompressible Brownian flow *<sup>g</sup>u*(*t*) with covariance *<sup>a</sup>* ∈ <sup>Γ</sup>(*T*O ⊗ *<sup>T</sup>*O) and time dependent drift *<sup>u</sup>*(*t*, ·) ∈ <sup>Γ</sup>(*T*O). We assume that for all *<sup>x</sup>* ∈ O, *<sup>a</sup>*(*x*, *<sup>x</sup>*) = <sup>2</sup>*ν***g**−1(*x*) for some *<sup>ν</sup>* > 0, where **<sup>g</sup>** denotes the metric tensor of the manifold. Such incompressible flows are known to be well defined on compact symmetric spaces and on compact Lie groups.

This means that

$$d\operatorname{g}\_{\mathfrak{u}}(t)(\mathfrak{x}) \otimes d\operatorname{g}\_{\mathfrak{u}}(t)(\mathfrak{y}) = a\left(\operatorname{g}\_{\mathfrak{u}}(t)(\mathfrak{x}), \operatorname{g}\_{\mathfrak{u}}(t)(\mathfrak{y})\right)dt,$$

$$d\operatorname{g}\_{\mathfrak{u}}(t)(\mathfrak{x}) \otimes d\operatorname{g}\_{\mathfrak{u}}(t)(\mathfrak{x}) = 2\nu \mathbf{g}^{-1}\left(\operatorname{g}\_{\mathfrak{u}}(t)(\mathfrak{x})\right)dt,$$

the drift of *gu*(*t*)(*x*) is absolutely continuous and satisfies *Dgu*(*t*)(*x*) = *u*(*t*, *gu*(*t*)(*x*)). Then, using the same kind of variations as in the flat case, we derive (c.f [23]), from the energy functional

$$S(\mathcal{g}^{\mu}) = \frac{1}{2}E\left[\int\_{0}^{T} \left(\int\_{\mathcal{O}} |D\mathcal{g}\_{\mu}(t)(\mathbf{x})|^{2} \, dm(\mathbf{x})\right) \, dt\right],$$

the equation

$$\frac{\partial}{\partial t}\mu + \nabla\_{\boldsymbol{\mu}}\mu = \nu \boldsymbol{L}\boldsymbol{\mu} - \nabla p$$

where *L* = *dd*∗ + *d*∗*d* is the the Laplace–de Rham operator. We recall that when computed on forms and, in particular, on vector fields, *L* differs from the usual Levi–Civita Laplacian by a Ricci curvature term.

#### *4.2. Lagrange Multipliers Approach*

Here we describe another stochastic variational principle for the Navier–Stokes equation, that uses a Lagrange multiplier formulation. Its advantage is that it does not need to be formulated in the space *GV*, since incompressibility is not incorporated in the definition of the flows, but given instead by the multiplier condition. The problem becomes strictly finite-dimensional and variations of the paths may be defined simply by shifts. It also allows us to consider domains with boundary. References for the Lagrange multipliers approach are [26,27], where one can find detailed proofs. The first concerns the case of a domain without boundary (more precisely the torus) and in the second one can find its extension to a domain with boundary.

Let O be a domain with a regular boundary. We consider the Navier–Stokes equation with Neumann boundary condition on [0, *T*] × O, namely the condition ∇*u*.*n* = 0, where *n* denotes the unit vector normal to the boundary. The stochastic Lagrangian flows can now take into account reflections at the boundary of the domain and can be written,

$$d\mathfrak{g}\_t(\mathbf{x}) = \sqrt{2\nu}d\mathcal{W}\_t + \mathfrak{u}(t, \mathfrak{g}\_t(\mathbf{x}))dt + \mathfrak{n}(\mathfrak{g}\_t(\mathbf{x}))d\ell(t), \qquad \mathfrak{g}\_0(\mathbf{x}) = \mathbf{x} \in \mathcal{O} \tag{8}$$

where (*t*) = *<sup>t</sup>* <sup>0</sup> 1*δ*O(*gs*(*x*))*d*(*s*) is the local time, representing the amount of time spent by the diffusion process in the neighborhood of points in the boundary.

The action functional is defined as

$$S(\mathbf{g}, p) = \frac{1}{2} \mathbb{E} \int\_0^T \int |D\_l \mathbf{g}\_l(\mathbf{x})|^2 dt dm(\mathbf{x}) + E \int\_0^T \int p(t, \mathbf{g}\_l(\mathbf{x})) (\det \nabla \mathbf{g}\_l(\mathbf{x}) - 1) dt dm(\mathbf{x}) \tag{9}$$
 
$$\coloneqq S^1(\mathbf{g}, p) + S^2(\mathbf{g}, p) \tag{10}$$

The extra term *S*<sup>2</sup> corresponds to a Lagrange multiplier whose constraint forces the paths to keep the volume measure preserved during the evolution (incompressibility condition). The variable *p* is defined in the linear space *<sup>L</sup>*2([0, *<sup>T</sup>*] × O). We consider variations of the form

$$\mathfrak{g}\_t(\cdot) \to \mathfrak{g}\_t^\epsilon(\cdot) = \mathfrak{g}\_t(\cdot) + \epsilon h(t, \mathfrak{g}\_t(\cdot))$$

$$p(t, \cdot) \to p^{\epsilon}(t, \cdot) = p(t, \cdot) + \epsilon \,\varphi(t, \gnot \cdot) \,\|$$

with *h*(*t*, *x*) and *ϕ*(*t*, *x*) deterministic and smooth, satisfying *h*(*T*, ·) = *h*(0, ·) = 0, *h* = 0 on *∂*O.

Concerning *S*2, we have,

$$\frac{d}{d\varepsilon}\Big|\_{\varepsilon=0}S^2(\mathbf{g}^\varepsilon, p^\varepsilon) = E \int\_0^T \int \varphi(t, \mathbf{g}\_t(\mathbf{x})) (\det \nabla \mathcal{g}\_t(\mathbf{x}) - 1) dt dm(\mathbf{x})\tag{11}$$

$$+E\int\_{0}^{T}\int\_{\underline{0}} (\nabla p(t,\underline{g}\_{t}(\mathbf{x})).h(t,\underline{g}\_{t}(\mathbf{x})) (\det \nabla \underline{g}\_{t}(\mathbf{x}) - 1)dt dm(\mathbf{x})\tag{12}$$

$$+E\int\_{0}^{T}\int p(t,\mathbf{g}\_{t}(\mathbf{x}))\,\frac{d}{d\varepsilon}\Big|\_{\varepsilon=0}\det\nabla(\mathbf{g}\_{t}(\mathbf{x})+\epsilon h(t,\mathbf{g}\_{t}(\mathbf{x}))dtdm(\mathbf{x})\tag{13}$$

Since *ϕ* is arbitrary we conclude from the first term of Equation (11) that critical points of the action are volume-preserving diffeomorphisms (det ∇*gt*(*x*) = 1) and therefore have divergence-free drifts. It follows immediately that Equation (12) is also equal to zero. The computation of the third term gives

$$\frac{d}{d\varepsilon}\Big|\_{\varepsilon=0}S^2(\mathfrak{g}^\varepsilon,p^\varepsilon)=-E\int\_0^T\int (\nabla p(t,\mathfrak{g}\_t(\mathbf{x}))\cdot h(t,\mathfrak{g}\_t(\mathbf{x})))dtdm(\mathbf{x}).$$

Actually this last computation does not use stochastic calculus:

$$\begin{aligned} \frac{d}{d\varepsilon}\Big|\_{\varepsilon=0} \det \nabla(\mathcal{g}\_t(\mathbf{x}) + \epsilon h(t, \mathcal{g}\_t(\mathbf{x}))) &= \det \nabla \mathcal{g}\_t(\mathbf{x}) \operatorname{tr} \Big( (\nabla \mathcal{g}\_t(\mathbf{x}))^{-1} \left. \frac{d}{d\varepsilon}\Big|\_{\varepsilon=0} \nabla \mathcal{g}\_t^{\varepsilon}(\mathbf{x}) \right) \\\\ &= \det(\nabla \mathcal{g}\_t(\mathbf{x})) \operatorname{tr} \Big( (\nabla \mathcal{g}\_t(\mathbf{x}))^{-1} \nabla (h(t, \mathcal{g}\_t(\mathbf{x}))) \Big). \end{aligned}$$

Since

$$
\partial\_i(p(t, \mathcal{g}\_t)(\nabla \mathcal{g}\_t)\_{ij}^{-1} h(t, \mathcal{g}\_t)^j) = \partial\_i(p(t, \mathcal{g}\_t)) (\nabla \mathcal{g}\_t)^{-1}\_{ij} h(t, \mathcal{g}\_t)^j + p(t, \mathcal{g}\_t) (\nabla \mathcal{g}\_t)^{-1}\_{ij} \partial\_i(h^j(t, \mathcal{g}\_t))
$$

$$
+ p(t, \mathcal{g}\_t) h^j(t, \mathcal{g}\_t) \partial\_i(\nabla \mathcal{g}\_t)^{-1}\_{ij},
$$

$$\mathbf{u}(13) = -E \int\_0^T \int [\partial\_i(p(t, \mathbf{g}\_t)) (\nabla \mathbf{g}\_t)^{-1}\_{ij} + p(t, \mathbf{g}\_t) \partial\_i((\nabla \mathbf{g}\_t)^{-1}\_{ij})] h^j(t, \mathbf{g}\_t) \det \nabla \mathbf{g}\_t \, dt d\mathbf{x}.$$

Notice that we already concluded that det ∇*gt* = 1. On the other hand,

$$\sum\_{i} \partial\_{i} (\nabla g\_{t})^{-1}\_{ij} = 0.$$

Indeed, derivating the equality det ∇*gt* = 1, we get

$$\partial\_k \det(\nabla \mathcal{g}\_t) = \text{tr}\left((\nabla \mathcal{g}\_t)^{-1} \partial\_k (\nabla \mathcal{g}\_t)\right) = \sum\_i (\nabla \mathcal{g}\_t)^{-1}\_{ij} \partial\_k \partial\_i \mathcal{g}^j\_t = 0.1$$

Also, derivating equality (∇*gt*)−<sup>1</sup> *ij ∂kg j <sup>t</sup>* = *δik*, we obtain

$$\sum\_{i} \partial\_{i} (\nabla \mathcal{g}\_{t})^{-1}\_{ij} \partial\_{k} \mathcal{g}^{j}\_{t} + (\nabla \mathcal{g}\_{t})^{-1}\_{ij} \partial\_{i} \partial\_{k} \mathcal{g}^{j}\_{t} = 0;$$

therefore

$$\sum\_{i} \partial\_{i} (\nabla \mathcal{g}\_{t})^{-1}\_{ik} = -\left( (\nabla \mathcal{g}\_{t})^{-1}\_{ij} \partial\_{k} \partial\_{i} \mathcal{g}^{j}\_{t} \right) (\nabla \mathcal{g}\_{t})^{-1}\_{jk} = 0$$

and

$$\mathbf{u}(13) = -\mathrm{E} \int\_0^T \int (\partial\_i(p(t, \mathcal{g}\_t(\mathbf{x}))) (\nabla \mathcal{g}\_t(\mathbf{x}))^{-1}\_{ij}) h^j\_t(\mathcal{g}\_t(\mathbf{x})) \det \nabla \mathcal{g}\_t(\mathbf{x})) dt dm(\mathbf{x})$$

*Water* **2020**, *12*, 864

$$=-E\int\_0^T \int (\nabla p(t, \mathcal{g}\_t(\mathbf{x})).h(t, \mathcal{g}\_t(\mathbf{x})))dt dm(\mathbf{x})\,.$$

We now look at the derivation of *S*1. We can prove, using Itô calculus and similarly to the computation in Section 4.1, that

$$\frac{d}{d\varepsilon}\Big|\_{\varepsilon=0}\mathcal{S}^{1}(\mathcal{g}^{\varepsilon},\mathcal{p}^{\varepsilon})=E-E-E\int\_{0}^{T}\int (D\_{t}D\_{t}g\_{t}(\mathbf{x})h(t,\mathcal{g}\_{t}(\mathbf{x})))dt d\mathbf{x}$$

$$-E\int\_{0}^{T}\int (dD\_{t}g\_{t}(\mathbf{x}).dh(t,\mathcal{g}\_{t}(\mathbf{x})))d\mathbf{x}.$$

$$=-E\int\_{0}^{T}\int (D\_{t}D\_{t}g\_{t}(\mathbf{x})h(t,\mathcal{g}\_{t}(\mathbf{x})))dt d\mathbf{m}(\mathbf{x})-2\nu E\int\_{0}^{T}(\nabla\mathbf{u}.\nabla h)(t,\mathcal{g}\_{t}(\mathbf{x}))dt d\mathbf{m}(\mathbf{x})$$

$$\text{With }\nabla h=\mathcal{H}\_{0},\mathcal{P}\_{0},h(\mathbf{x})=(\partial\_{t}+\nu(\mathbf{x}))\text{ and }\nabla\mathcal{H}\_{0},\mathcal{P}\_{0},h(\mathbf{x})=\partial\_{t}\mathcal{H}\_{0},\text{ we have }$$

Using equality *DtDtgt*(*x*)=( *<sup>∂</sup> <sup>∂</sup><sup>t</sup> u* + (*u*.∇)*u* + *ν*Δ*u* + *n*∇*u* (*t*))(*t*, *gt*(*x*)) and integration by parts, we deduce that

$$\frac{d}{d\varepsilon}\Big|\_{\varepsilon=0} \mathcal{S}^1(\mathcal{g}^\varepsilon, p^\varepsilon) = -E \int\_0^T \int [((\frac{\partial}{\partial t}u + (u.\nabla)u - \nu \Delta u), h) + (n.\nabla u) \ h \,\ell(t)](t, \mathcal{g}\_t(x))dt dm(x).$$

Combining the expressions above for the variation of the action functional and using the invariance of the volume measure for the flows, we obtain the following result,

**Theorem 2.** *A diffusion gt of the form of Equation (8) and a function p are critical for the action functional in Equation (9) if and only if the drift u*(*t*, ·) *of gt satisfies the Navier–Stokes equation*

$$
\partial\_t u + (u.\nabla)u = \nu \Delta u - \nabla p, \qquad \operatorname{div} u(t, \cdot) = 0, \quad \nabla u.\mathbb{n} = 0 \text{ in } \partial \mathcal{O} \tag{14}
$$

*with t* ∈ [0, *T*]*.*

#### **5. Constructing Solutions of Navier–Stokes Equation by Probabilistic Methods**

#### *5.1. Forward–Backward Stochastic Differential Systems*

The stochastic Lagrangian flows *gt*(·) obtained either via the geometric or the Lagrangian multiplier approach described above satisfy an equation of Newton type. More precisely, on a fixed time interval [0, *T*] and after the change of time *v*(*t*, *x*) = −*u*(*T* − *t*, *x*), we have *Dgt*(*x*) = −*v*(*T* − *t*, *gt*(*x*)) and

$$D\_t v(t, \mathfrak{g}\_t(\mathbf{x})) = -\nabla p(t, \mathfrak{g}\_t(\mathbf{x})), \qquad \mathfrak{g}\_0(\mathbf{x}) = \mathbf{x} \tag{15}$$

or, regarded as an equation in the weak *L*<sup>2</sup> sense,

$$D\_l v(t, \mathfrak{g}\_l(\cdot)) = 0, \qquad \mathfrak{g}\_0 = \mathrm{id}.\tag{16}$$

Formally, when the viscosity is zero, this is Arnold's geodesic equation and ours can be seen, indeed, as a generalized geodesic equation. Moreover, although we only consider in this paper an Euclidean setting, the framework can be extended without essential difficulties to a general Riemannian manifold or a Lie group ([23,24]).

How does one solve directly this equation of motion? One possible way is to characterize it in terms of forward–backward stochastic differential systems, which are second order stochastic equations (c.f the Appendix A). We have formulated the problem in [28] first, then solved it in [29,30] for some specific function spaces and in two and three dimensions, respectively. The characterization in terms of forward–backward stochastic differential equations has also the advantage that it may allow to implement numerical methods, known for such systems (c.f for example, [31] and the recent work [32]).

The forward–backward system solved by our stochastic Lagrangian flows can be written in the form

$$\begin{cases} \,d\mathcal{g}\_t = \sqrt{2\nu} \,dW\_t + \mathcal{Y}\_t dt \\ dY\_t = Z\_t dW\_t - \nabla p(\mathcal{g}\_t) dt \end{cases}$$

together with a given initial condition for the forward equation (*g*<sup>0</sup> = *x*) and a final condition for the backward one, *YT*. It can be proved these type of systems are well posed and their solutions are of the form *Yt* = *Yt*(*x*) = *v*(*t*, *gt*(*x*))(= −*u*(*t*, *gt*(*x*))) for some vector field *v*. We have,

$$D\_t \mathbf{\hat{y}}\_t = D\_t D\_t \mathbf{g}\_t = -\nabla p(\mathbf{g}\_t)$$

(compare with Equation (15)). Variable *Z*, although a priori unknown, is a posteriori determined and equal to <sup>√</sup>2*ν*∇*v*(*gt*).

To be more precise, let *u* be a solution of the Navier–Stokes equation in the time interval *<sup>t</sup>* ∈ [0, *<sup>T</sup>*] and assume that *<sup>u</sup>* is regular. Let *<sup>g</sup><sup>t</sup> <sup>s</sup>*(*x*) be the unique solution of the following stochastic differential equation,

$$\begin{cases} \quad d\mathfrak{g}\_s^t(\mathbf{x}) = \sqrt{2\nu} d\mathcal{W}\_s - \mathfrak{u}(T - s\_\prime \mathfrak{g}\_s^t(\mathbf{x}))ds \\\quad \mathfrak{g}\_t^t(\mathbf{x}) = \mathfrak{x}, \end{cases}$$

with *s* > *t*. We define *Y<sup>t</sup> <sup>s</sup>*(*x*) = *<sup>u</sup>*(*<sup>T</sup>* − *<sup>s</sup>*, *<sup>g</sup><sup>t</sup> <sup>s</sup>*(*x*)), *Z<sup>t</sup> <sup>s</sup>*(*x*) = ∇*u*(*<sup>T</sup>* − *<sup>s</sup>*, *<sup>g</sup><sup>t</sup> <sup>s</sup>*(*x*)); applying Itô's formula directly, the following forward–backward stochastic differential system with solution (*g<sup>t</sup> s*(*x*),*Y<sup>t</sup> <sup>s</sup>*(*x*), *Z<sup>t</sup> <sup>s</sup>*(*x*), *u*(*t*, *x*), *p*(*t*, *x*)) is derived,

$$\begin{cases} & d\mathcal{g}\_s^t(\mathbf{x}) = \sqrt{2\nu}dW\_s - u(T - s, \mathcal{g}\_s^t(\mathbf{x}))ds \\ & dY\_s^t(\mathbf{x}) = \sqrt{2\nu}Z\_s^t(\mathbf{x})dW\_s - \nabla p(T - s, \mathcal{g}\_s^t(\mathbf{x}))ds \\ & Y\_t^t(\mathbf{x}) = u(T - t, \mathbf{x}) \\ & \mathcal{g}\_t^t(\mathbf{x}) = \mathbf{x}, \ Y\_T^t(\mathbf{x}) = u\_0(\mathcal{g}\_T^t(\mathbf{x})) \end{cases} \tag{17}$$

Recall also that the pressure satisfies the following identity

$$
\Delta p(t, \mathbf{x}) = \sum\_{i, j = 1}^{3} \partial\_i u^j(t, \mathbf{x}) \partial\_j u^i(t, \mathbf{x}). \tag{18}
$$

On the other hand, if (*g<sup>t</sup> s*(*x*),*Y<sup>t</sup> <sup>s</sup>*(*x*), *Z<sup>t</sup> <sup>s</sup>*(*x*), *u*(*t*, *x*), *p*(*t*, *x*)) is a solution of Equation (17) together with Equation (18) and *<sup>u</sup>* is regular enough, then the vector field *<sup>u</sup>*(*t*, *<sup>x</sup>*) :<sup>=</sup> *<sup>Y</sup>T*−*<sup>t</sup> <sup>T</sup>*−*t*(*x*) satisfies the Navier–Stokes equation for *t* ∈ [0, *T*]. In particular, we can show that div *u*(*t*, *x*) = 0 due to the expression of *p*(*t*, *x*) given by Equation (18).

In order to incorporate Equation (18) in the forward–backward system and obtain a closed system of equations we proceed as follows. Denote by *<sup>N</sup>* the Newton's potential in <sup>R</sup>*d*, *<sup>d</sup>* <sup>≥</sup> 3, i.e., the operator Δ−<sup>1</sup> which is given by

$$Nf(\mathbf{x}) = \mathcal{C}(d) \int\_{\mathbb{R}^d} \frac{f(y)}{|\mathbf{x} - y|^{d-2}} dy \,\mathrm{s}$$

where *C*(*d*) is a constant depending on the dimension *d*. Then we can write Equations (17) and (18) as

$$\begin{cases} & d\mathcal{g}\_s^t(\mathbf{x}) = \sqrt{2\nu} d\mathcal{W}\_s - u(T - s, \mathcal{g}\_s^t(\mathbf{x})) ds \\ & dY\_s^t(\mathbf{x}) = \sqrt{2\nu} Z\_s^t(\mathbf{x}) d\mathcal{W}\_s - \nabla N\Big(\sum\_{i,j} \partial\_i \boldsymbol{\upsilon}^j - \partial\_j \boldsymbol{\upsilon}^j\Big) (T - s, \mathcal{g}\_s^t(\mathbf{x})) ds \\ & Y\_t^t(\mathbf{x}) = u(T - t, \mathbf{x}) \\ & \mathcal{g}\_t^t(\mathbf{x}) = \mathbf{x}, Y\_T^t(\mathbf{x}) = \mathbf{u}\_0(\mathcal{g}\_T^t(\mathbf{x})), \end{cases} \tag{19}$$

Using suitable *<sup>L</sup><sup>p</sup>* bounds of the operator ∇*<sup>N</sup>* <sup>∑</sup>*i*,*<sup>j</sup> <sup>∂</sup>iv<sup>j</sup>* − *<sup>∂</sup>jv<sup>i</sup>* we have constructed in [30] local unique solutions of the system in Equation (19), and therefore of the Navier–Stokes equation, in some Sobolev-type functional spaces in dimension *d* ≥ 3. The two-dimensional case was studied in a similar way [29], but via the vorticity equation.

Such type of forward–backward differential equations were also studied on general Lie groups in [33].

#### *5.2. Entropy Methods*

In a closely related recent work ([34]), we explored the fact that the stochastic kinetic energy used here coincides in fact with a relative entropy, of the law of the diffusion process associated with Navier–Stokes in relation to the Wiener law. Brenier's relaxation of Arnold's geodesic problem, extended to the stochastic setting, becomes, in this sense, an entropy minimization problem.

Recall that the relative entropy of a probability measure *Q* on some measurable set *X* with respect to a reference probability measure *R* on *X* is defined as

$$H(Q|R) = \int\_{\mathcal{X}} \log\left(\frac{dQ}{dR}\right) dQ$$

when *Q* is absolutely continuous with respect to *R*. What we have named the Brenier–Schrödinger problem in [34] consists of the following. Consider *X* to be an Euclidean space and fix as reference measure *R* as the law of the reversible Brownian motion on a time interval [0, *T*], namely *R* = *<sup>X</sup> <sup>P</sup>xdm*(*x*) for *<sup>P</sup><sup>x</sup>* the law of a Brownian motion starting at *<sup>x</sup>*. The we look for a measure *<sup>Q</sup>* minimizing the relative entropy *H*(*Q*|*R*) under the constraint that *Qt*, the marginals at time *t*, are given, and that the joint law at initial and final times, *Q*0*T*, is also given. In our case the relevant condition for the marginals is the incompressibility condition *Qt* = *m* for all times.

This is a convex minimization problem reminiscent of the Schrödinger problem in Quantum Mechanics ([35–37]) and, simultaneously, a generalization of a Monge–Kantarovich mass transportation problem (c.f [38]).

The relation between the entropy problem and our stochastic variational one relies on the fact that, by Girsanov's theorem (c.f for instance [19]), if the Lagrangian path is given by the stochastic differential equation *dg<sup>u</sup> <sup>t</sup>* <sup>=</sup> <sup>√</sup>2*νdWt* <sup>+</sup> *<sup>u</sup>*(*t*, *<sup>g</sup><sup>u</sup> <sup>t</sup>* )*dt* with initial condition *R*, and if *Q* denotes the law of this stochastic process, we have

$$H(\mathbf{Q}|\mathbf{R}) = H(\mathbf{Q}\_0|\mathbf{R}) + \frac{1}{2}E \int\_0^T |u\left(t, \mathbf{g}\_t^u\right)|^2 dt$$

Therefore the entropy minimization problem corresponds to our stochastic variational principle, described here in Section 4. We can use then tools of convex analysis, establish the existence of a solution, consider a dual problem that provides the pressure term and exhibit Lagrange multipliers. Since our solution is the law of a stochastic process which is absolutely continuous with respect to the reference measure, one can write its Radon–Nikodym derivative and compare it to an alternative expression obtained from convex duality. We find that the drifts of our Lagrangian flows, written as velocity fields, have unusual properties: when stratified (conditioned to the arrival point of the stochastic process), the backward velocity is the gradient of the solution to a Hamilton–Jacobi–Bellman equation, while the pressure remains independent of the arrival point. We recover, with this method, a structure that is analogous to Brenier's stratified solutions in [39].

#### *5.3. A Weak Notion of Navier–Stokes Solutions*

One cannot expect that the stochastic variational approach will determine solutions of the Navier–Stokes in full generality. It is already the case that geodesics associated with Euler equation do not exist in some situations. The main difficulty is that the topology induced by the energy is not strong enough to deal with the needed regularity of the Lagrangian maps. In order to overcome the problem, Brenier has introduced a concept that he named "generalized solutions", replacing the notion of geodesic path by a probability measure over geodesic paths, in the spirit of the Monge–Kantarovich problem ([40]). An extension of Brenier's framework to the stochastic Lagrangian setting (and in particular to the Navier–Stokes equation) can be found in [41].

For the Euler equation, Brenier looked for probability measures *P* on the path space *C*([0, *T*]; O) which minimize the energy functional

$$\frac{1}{2} \int\_{C([0,T];\mathcal{O})} \left[ \int\_0^T |\dot{\gamma}(t)|^2 dt \right] dP(\gamma)$$

with the incompressibility constraints (*et*)∗*P* = *dm*(*x*), where *et* : *γ* → *γ*(*t*) is the evaluation map. A solution of such problem gives rise to a weaker type of Lagrangian flows for the Euler equation. Indeed one can define a probability measure *<sup>μ</sup>* on [0, *<sup>T</sup>*] × O<sup>2</sup> by

$$\int f(t,x,v)d\mu = \frac{1}{T} \int\_{C([0,T];\mathcal{O})} f(t,\gamma(t),\gamma'(t))dP(\gamma)dt$$

and *μ* solves the Euler equation in a weak sense, namely

$$\int [\upsilon . w'(x)a'(t) + \upsilon . (\nabla w(x). \upsilon)a(t)] d\mu(t, x, \upsilon) = 0$$

holding for all smooth test functions *α*(*t*) and every smooth divergence free vector field *v*. These kind of weak solutions of partial differential equations is known as solutions in the sense of Di Perna and Majda.

In the viscous case, we consider <sup>O</sup>-valued semimartingales *gt* of the form *dgt* <sup>=</sup> <sup>√</sup>2*νdWt* <sup>+</sup> *utdt* (remark that the drift can be random) and their corresponding laws <sup>P</sup>*<sup>g</sup>* on the path space *<sup>C</sup>*([0, *<sup>T</sup>*]; <sup>O</sup>): for every cylindrical functional *F*,

$$\int\_{C([0,T],\mathcal{O})} F(\gamma(t\_1), \dots, \gamma(t\_n)) d\mathbb{P}^{\mathbb{g}}(\gamma) = \int\_{\mathcal{O}} \left[ \int\_{C([0,T],\mathcal{O})} F(g\_{t\_1}(\mathbf{x}), \dots, g\_{t\_n}(\mathbf{x})) d\mathbb{P}^{\mathbb{g}}\_{\mathbf{x}} \right] dm(\mathbf{x}) \ge 0$$

where P*<sup>g</sup>* = P*<sup>g</sup> <sup>x</sup>* <sup>⊗</sup> *dm*(*x*). Under <sup>P</sup>*<sup>g</sup> <sup>x</sup>*, the semimartingale *gt* starts from *x*.

We say that the semimartingale *gt* is incompressible if, for each *t* > 0,

$$E\_{\mathbb{P}^\sharp}[f(\mathcal{g}\_t)] = \int\_{\mathcal{O}} f(\mathbf{x})d\mathbf{x}, \quad \text{for all } f \in \mathcal{C}(\mathcal{O}).$$

and we define the energy functional of a semimartingale *g* by the formula

$$\frac{1}{2} \int\_{\mathcal{O}} E\_{\mathbb{P}^{\mathbb{E}}\_{x}} (\int\_{0}^{T} |D\_{t}g(\mathbf{x})|^{2} dt) dm(\mathbf{x}) = \frac{1}{2} \int\_{\mathcal{O}} E\_{\mathbb{P}^{\mathbb{E}}\_{x}} \left( \int\_{0}^{T} |u\_{t}(g\_{t}(\mathbf{x}))|^{2} dt \right) dm(\mathbf{x})$$

Considering suitably admissible variations we concluded in [41] that the semimartingales which are critical points of the energy functional above solve the Navier–Stokes equation in the following weak form:

$$\int\_{\mathcal{O}} \int\_{0}^{T} < \boldsymbol{u}\_{l} \, a'(t)\boldsymbol{w} + a(t) \, \nabla \boldsymbol{w} \cdot \boldsymbol{u}\_{l} - \nu \, a(t) \Delta \boldsymbol{w} > \, \mathrm{d}t \mathrm{d}m(\mathbf{x}) = 0 \tag{20}$$

for all smooth functions of time *α* and all smooth vector fields *w* such that div(*w*) = 0. Moreover we have showed that, in the case where O is the torus (corresponding to periodic boundary conditions) and under certain additional assumptions, classical solutions of the Navier–Stokes are minimizers of the energy action functional.

#### **6. Stability Properties**

In finite dimensions it is well known that the behavior of geodesics can be expressed in terms of the curvature of the underlying manifold via the Jacobi equation.

Arnold's approach to the Euler equation allowed to show, in many cases, that the curvature of the spaces of diffeomorphisms is negative and therefore that the fluid trajectories are unstable (or "chaotic"), i.e., their distance, starting from different initial conditions, grows exponentially during time evolution (c.f [16]).

For viscous flows it is expected that particles become closer and closer after some possible initial stretching. For our model, at least in the case of the two dimensional torus, we could show ([23] and also [42]) that sensitivity with respect to initial conditions of the trajectories is enhanced by their stochasticity. The behavior will depend of the choice of diffusion coefficients that we consider in the Lagrangian paths or, in other words, on which scales and with what strength the motion is excited.

We have considered on the two-dimensional flat torus T = (R/2*π*Z)<sup>2</sup> a Brownian motion of the form

$$dB\_t = \sum\_{k \in \mathbb{Z}^2} (k\_{2\prime} - k\_1) \sqrt{\nu} \lambda\_k A\_k dW\_t^k$$

where *<sup>A</sup>*2*k*(*x*) = cos(2*k*.*x*), *<sup>A</sup>*2*k*+1(*x*) = sin((2*<sup>k</sup>* <sup>+</sup> <sup>1</sup>).*x*), *<sup>k</sup>* <sup>∈</sup> <sup>Z</sup>2, *<sup>x</sup>* <sup>∈</sup> <sup>T</sup>. The coefficients *<sup>λ</sup><sup>k</sup>* have to be such that the process *B* is well defined, which is achieved by assuming a certain decay at infinity of those coefficients (see [22]). The Lagrangian stochastic processes associated to the Navier–Stokes equation will be defined as in Equation (7), with *Wt* replaced by *Bt*.

Consider the following distance between trajectories:

$$
\rho^2(\mathcal{g}, \tilde{\mathfrak{g}}) = \int\_{\mathbb{T}} |\mathfrak{g}(\mathfrak{x}) - \tilde{\mathfrak{g}}(\mathfrak{x})|^2 dm(\mathfrak{x}).
$$

By using Itô calculus we have deduced, in particular, that under the assumption that the Navier–Stokes solution *<sup>u</sup>* satisfies ∇*u*(*t*, *<sup>x</sup>*) ≤ *<sup>c</sup>*1*e*−*c*2*<sup>t</sup>* , we have an estimate for the distance between trajectories of the form

$$
\rho\_t \ge \rho\_5 \exp\left(Z\_t + c\_3 t - \frac{c\_1}{c\_2} (1 - e^{-c\_2 t})\right)
$$

for *s* ≤ *t*, where *Z* is a 1-dimensional Brownian motion, *c*<sup>3</sup> another constant which depends on the coefficients *λ*.

The assumption on the gradient of *u* implies that the velocity decays to zero at exponential rate. On the contrary the stochastic Lagrangian flows, describing the position of the fluid, get apart exponentially, at least for short times. Moreover, by the explicit expression of the constant *c*<sup>3</sup> ([23]), we observe that the stochastic Lagrangian trajectories for a fluid with a given viscosity constant tend to get apart faster when the higher Fourier modes (and therefore the smaller length scales) are randomly excited.

We can also show how the rotation of two particles, when their distance is small, becomes more and more irregular as time evolves, with explicit formulae in the torus case.

#### **7. A Stochastic Navier–Stokes Equation**

When we consider random forces, Navier–Stokes becomes a stochastic partial differential equation. In this paragraph we formulate some of these type of equations by a variational principle, where the action functional is suitably perturbed.

*Water* **2020**, *12*, 864

We define the action functional, which is now a random variable (in particular we remove the expectation) with some extra stochastic terms in the Lagrangian, as

$$S\_{\omega}(\mathbf{g}, p) = \frac{1}{2} \int\_{0}^{T} \int |D\_{l}\mathbf{g}\_{l}(\mathbf{x})|^{2} dt dm(\mathbf{x}) + \int\_{0}^{T} \int p(t, \mathbf{g}\_{l}(\mathbf{x})) (\det \nabla \mathbf{g}\_{l}(\mathbf{x}) - 1) dt dm(\mathbf{x}) \tag{21}$$

$$+\int\_{0}^{T}\int D\_{l}\mathbf{g}\_{l}(\mathbf{x})d\mathcal{M}\_{l}^{\mathbb{S}}(\mathbf{x})dm(\mathbf{x}) - \sqrt{2\nu}\int\_{0}^{T}\int D\_{l}\mathbf{g}\_{l}(\mathbf{x})d\mathcal{W}\_{l}dm(\mathbf{x})\tag{22}$$

where *M<sup>g</sup>* denotes the martingale part of the diffusion process *g*. We use the same variations as in Section 4.2. The variation of the two extra terms gives

$$\int\_{0}^{T} \int [(h(t, \mathfrak{g}\_t), \sqrt{2\nu} d\mathcal{W}\_t) + (D\_t \mathfrak{g}\_t.(\nabla h(t, \mathfrak{g}\_t).d\mathcal{W}\_t)) - (h(t, \mathfrak{g}\_t).\sqrt{2\nu} d\mathcal{W}\_t)] dm(\mathbf{x}) = 0$$

that reduces to

$$\int\_0^T \int v(t, \mathbf{g}\_t(\mathbf{x}). (\nabla h(t, \mathbf{g}\_t(\mathbf{x})).dW\_t)d\mathbf{x} = \int\_0^T v(t, \mathbf{x}). (\nabla h(t, \mathbf{x}).dW\_t)$$

$$= -\int\_0^T ((\nabla v(t, \mathbf{x}).h(t, \mathbf{x})).dW\_t)$$

equality which holds for all *h*, *P*-almost surely.

Hence we obtain the following (c.f [25,26]),

**Theorem 3.** *A diffusion gt of Equation (8) and a function p are critical for the action functional Sω if and only if the drift u*(*t*, ·) *of gt satisfies the stochastic Navier–Stokes equation*

$$d\_l u + (u.\nabla)u dt = \sqrt{2\nu} \nabla u d\mathcal{W}\_l + \nu \Delta u dt - \nabla p dt, \qquad \text{div } u(t, \cdot) = 0, \quad \nabla u.n = 0 \text{ in } \partial \mathcal{O} \tag{23}$$

*with t* ∈ [0, *T*]*.*

The stochastic Navier–Stokes equation written above can also be regarded as a stochastic perturbation of Euler equation in the sense that, replacing Itô differentials by Stratonovich ones, Equation (23) is equivalent to

$$d\_t u + (u. \nabla) u dt = \sqrt{2\nu} \nabla u \circ d\mathcal{W}\_t - \nabla p dt, \qquad \text{div } u(t, \cdot) = 0, \quad \nabla u. n = 0 \text{ in } \partial \mathcal{O}$$

with *t* ∈ [0, *T*] and where ◦*d* stands for Stratonovich differential.

We remark that much more general noises can be considered and that we have only chosen a simple Brownian motion for simplicity. On the other hand the fact that we obtain a transport type noise is intrinsically related to our model.

There are many references for stochastic partial differential equations which are perturbations of Euler or Navier–Stokes. One can find in the recent work [5] a model where, as in our case, the noise is multiplicative.

In a general Lie group framework a variational approach to stochastic partial differential equations was developed in [25].

#### **8. Other Equations, Methods and Discussion**

The stochastic approach to Navier–Stokes described in this review can, in principle, be applied to describe any second order perturbation of a conservative first-order ordinary differential equation, when the latter is formulated in a suitable geometric way, as the Euler equation is formulated via a

geodesic equation. This is the case, for example, for the Camassa–Holm equation, the Hunter–Saxton equation, the average Euler equation and the equations governing the motion of rigid bodies.

The viscous Camassa–Holm equation was studied in detail in [43]. It corresponds to replacing the *L*2(*dm*) metric, the energy in the Lagrangian, by the Sobolev *H*<sup>1</sup> metric and this equation reads

$$\frac{\partial}{\partial t}v + (\mu . \nabla)v = \nu \Delta v + \sum\_{j} \nabla \mu^j. \Delta \mu^j - \nabla p\_{\prime} \quad \text{div} \boldsymbol{\mu} = 0$$

where *v* = *u* − Δ*u*.

Compressible Navier–Stokes or viscous Magnetohydrodynamical equations, on the other hand, have to be coupled with tracers (advected quantities) that need to be modeled by extra variables. Mathematically one introduces, together with the Lie group for the Lagrangian flows, a semidirect product vectorial structure (c.f [44]). We refer to [25] for a study of these various equations where we use stochastic methods, and stress again the fact that we can derive such equations by means of variational principles without any reference to thermodynamic considerations.

An extension of Arnold's geometric framework, describing compressible fluid dynamics (and other systems) with Newton's equations on a space of probability densities, can be found in the recent paper [45]. It reveals interesting connections with optimal transport and information theory, in particular.

The construction of solutions as critical points of the stochastic action functionals does not easily follow from the direct methods of the calculus of variations. We have instead indicated an indirect approach based on forward–backward stochastic differential equations, but even there we encountered many technical difficulties. Other possibilities are the use of entropy methods, as described in Section 5.2, or the relaxation of the notion of solution (c.f Section 5.3).

For conservative dynamical systems, it is well known that one can take advantage of the symmetries of such systems to reduce the complexity of the equations. Symmetries also play an important role in the implementation of numerical algorithms designed to investigate the equations in question. The natural objects to be introduced in order to replace constants of motion are martingales, since a martingale *M*, by definition, is a quantity whose generalized derivative vanishes (*DtM* = 0). One can find a Noether theorem for a certain class of diffusion processes in [46]. We refer to [47] for a brief discussion on symmetries in our context.

Finally, let us point out that numerical methods tied up with our probabilistic approach still have to be developed.

**Funding:** This research was partly funded by FCT (Fundação para a Ciência e Tecnologia, Portugal), grant "Schrödinger's problem and Optimal Transport: a multidisciplinary perspective", with reference PTDC/MAT-STA/28812/2017.

**Conflicts of Interest:** The author declare no conflict of interest.

#### **Appendix A. Some Basic Notions of Stochastic Calculus**

We start by recalling the definition of a (R-valued) *Brownian motion Wt*, *<sup>t</sup>* <sup>∈</sup> [0, *<sup>T</sup>*]. This is a continuous stochastic process (namely, a function of time and of a random parameter in some standard probability space (Ω, B, *P*)) that verifies the following properties:

(i) *W*<sup>0</sup> = *x*;


A R*<sup>d</sup>* Brownian motion corresponds to a collection of independent and identically distributed copies of one dimensional Brownian motions. A Brownian motion is also called a Wiener process.

We generally omit the random variable in the notation when we write a stochastic process.

Suppose that the probability space is endowed with a filtration, namely an increasing family P*<sup>t</sup>* of sub *σ*-algebras of B. Typically each P*<sup>t</sup>* represents the events that occur before a time *t*. A stochastic process *Mt* is a (R-valued) *martingale* with respect to <sup>P</sup>*<sup>t</sup>* if


A real-valued Brownian motion can also be characterized as a martingale with continuous sample paths such that, for all *t*, *W*<sup>2</sup> *<sup>t</sup>* − *t* is also a martingale.

When a martingale is (a.s.) continuous in time and satisfies the assumption *E*|*Mt*| <sup>2</sup> < +∞, we define its quadratic variation as the limit, in probability, of the sums

$$\_t = \sum\_{t\_i, t\_{i+1}} (M\_{t\_{i+1}} - M\_{t\_i})^2$$

when the mesh of the partition {*ti*} goes to zero. For a Brownian motion defined in [0, *T*] and starting from zero, this limit is equal to *T*. One defines the covariation between two martingales *Mt* and *Nt* as the limit of the sums

$$<\langle M, N>\_t = \sum\_{t\_i, t\_{i+1}} (M\_{t\_{i+1}} - M\_{t\_i})(N\_{t\_{i+1}} - N\_{t\_i})\dots$$

A stochastic process *Xt* is a *semimartingale* if, for every *t* it can be decomposed into a sum

$$\mathbf{X}\_t = M\_t + A\_t$$

of a martingale *Mt* and a process of bounded variation (or almost-surely differentiable) *At*, with *A*<sup>0</sup> = 0 a.s. The martingale part describes a diffusion, the bounded variation part an evolution which is "similar" to a deterministic one. The quadratic variation of a semimartingale is, by definition, equal to the quadratic variation of its martingale part.

Martingales and, in particular, Brownian motion, are (almost surely) never differentiable in time; therefore one cannot integrate with respect to them with a Lebesgue–Stieltjes kind of integration. Itô introduced a theory of stochastic integration, now named *Itô's stochastic calculus*. The Itô's integral is defined by the following limit

$$\int\_0^t X(s)dW(s) = \lim\_{t\_{i^\*}t\_{i+1}} X(t\_i) (\mathcal{W}\_{t\_{i+1}} - \mathcal{W}\_{t\_i})\_\*$$

the limit being taken in probability and when we consider partitions of the time interval [0, *t*] with mesh converging to zero. Itô's integral is well defined when *Xt* is a semimartingale such that *E*|*Xt*| <sup>2</sup>*dt* < +∞.

Notice that the values of the integrand *Xt* are taken on the left point of the intervals [*ti*, *ti*+1]. Unlike in usual Lebesgue–Stieltjes integration, considering its values in any other point of these intervals leads to completely different results. Another common and interesting way to define a stochastic integral is the so-called *Stratonovich integral*:

$$\int X(s) \diamond dW(s) = \lim\_{\mathcal{L}} \frac{1}{2} \sum\_{t\_i, t\_{i+1}} (X(t\_i) + X\_{t\_{i+1}}) (W\_{t\_{i+1}} - W\_{t\_i}).$$

Each type of integrals present its own advantages. Although Stratonovich integral demands more regularity on the integrand in order to be well defined, it is a more intrinsic concept and more adapted to be extended to semimartingales with values in curved spaces, for example. Also, as we shall write below, the rules of differential calculus are analogous to the classical ones for these integrals (and quite different for the Itô integral), which makes Stratonovich somehow popular in applications or in Physics. Nevertheless it turns out that *<sup>t</sup>* <sup>0</sup> *<sup>X</sup>*(*s*)*dW*(*s*) is a martingale, where *<sup>t</sup>* <sup>0</sup> *X*(*s*) ◦ *dW*(*s*) is not, in general. More important for us, any martingale is in fact a stochastic Itô integral and therefore any semimartingale can be decomposed into an Itô integral representing a diffusion, or a pure fluctuation, and a bounded variation part that represents the mean dynamical content of the process. So even though, when both integrals are defined, they are related by the following formula,

$$\int\_0^t X(s) \circ dW(s) = \int\_0^t X(s)dW(s) + \frac{1}{2} \int\_0^t d < X, W >\_{s,t}$$

dynamics is better identified using Itô integration.

Rules of differentiation are given by the so-called *Itô's formula*. If *f* is a regular function and *Xt* a (R*<sup>d</sup>* valued) semimartingale,

$$f(\mathbf{X}\_t) = f(\mathbf{X}\_0) + \sum\_{i=1}^d \int\_0^t \frac{\partial}{\partial \mathbf{x}\_i}(\mathbf{X}\_s) dX^i(\mathbf{s}) + \frac{1}{2} \sum\_{i,j=1}^d \int\_0^t \frac{\partial^2}{\partial \mathbf{x}\_i \partial \mathbf{x}\_j} f(\mathbf{X}\_s) d\mathbf{s} < \mathbf{X}^i, \mathbf{X}^j >\_{\mathbf{s}} \mathbf{0}$$

or, in differential form,

$$df(X\_t) = \sum\_{i=1}^d \frac{\partial}{\partial \mathbf{x}\_i}(X\_t)dX^i(t) + \frac{1}{2} \sum\_{i,j=1}^d \frac{\partial^2}{\partial \mathbf{x}\_i \partial \mathbf{x}\_j} f(X\_t)d < X^i, X^j > t \dots$$

For Stratonovich integrals we have, as in usual differential calculus,

$$df(X\_t) = \sum\_{i=1}^d \frac{\partial}{\partial \mathfrak{x}\_i}(X\_t) \circ dX^i(t).$$

If *f* is time dependent one has to add the term *<sup>∂</sup> <sup>f</sup> <sup>∂</sup><sup>t</sup>* (*t*, *Xt*)*dt* in the formulae.

Recall that, in the case of independent Brownian motions *W<sup>i</sup> <sup>t</sup>*, and writing id(*t*) = *t*, we have

$$d < \mathcal{W}^t, \mathcal{W}^l >\_t = \delta\_{ij} \, dt, \qquad d < \mathcal{W}^t, \text{id} >\_t = 0 \, \, \forall i.$$

*Stochastic differential equations* generalize ordinary differential ones. Given *σ* with values in matrices of type *d* × *r*, a vector field (possibly time-dependent) *u* and an initial random variable *X*0, these equations take the form

$$dX\_t = \sigma(X\_t) \, dW(t) + \mu(t, X\_t) dt.$$

Using Itô's formula we deduce that, if *Xt* is a solution of such an equation and *X*<sup>0</sup> = *x*, we have, for a smooth function *f* ,

$$\lim\_{t \to 0} \frac{1}{t} E\_x(f(X\_t) - f(0, x)) = \mathcal{L}f(x)\_r$$

where L is the second order linear operator

$$\mathcal{L}f(\boldsymbol{x}) = \frac{1}{2} \sum\_{i,j} (\sigma \boldsymbol{\sigma}^T)^{i}\_{i} \frac{\partial^2 f}{\partial \boldsymbol{\pi}\_i \partial \boldsymbol{\pi}\_j} + (\boldsymbol{\mu}, \boldsymbol{\nabla}f).$$

The operator L is called the generator of the process *Xt*.

Solutions of stochastic differential equations are particular cases of semimartingales. They are defined globally in time when the coefficients *σ* and *u* are smooth and bounded. Many works, even very recent ones, have been devoted to extend existence of solutions for less regular coefficients, in particular of Sobolev type.

When considering diffusion processes on domains rather than in all the Euclidean space, and subject to boundary conditions, we need to use the notion of *local time*. The local time of a Brownian motion is a family of (a.s.) continuous non negative random variables (*t*, *x*) such that, for any set *A* and *t* > 0 we have

$$\int\_0^t \mathbf{1}\_A(\mathcal{W}\_s) ds = 2 \int\_A \ell(t, x) dm(x).$$

It may also be regarded as the limit

$$\ell(t, x) = \lim\_{\varepsilon \to 0} \frac{1}{4\varepsilon} \int\_0^t \mathbf{1}\_{]\_{\varkappa - \varepsilon, x + \varepsilon[}(\mathcal{W}\_s) \, ds} .$$

Local time corresponds to the amount of time spent by a Brownian path in a neighborhood of a point *<sup>x</sup>* <sup>∈</sup> <sup>R</sup>. This concept was introduced by Paul Lévy in 1948 (see, for example [48]). The concept extends naturally to the multidimensional case as well as to semimartingales *Xt*. We have,

$$\int\_0^t f(X\_t)dt < X\_\prime X >\_t = 2\int\_{-\infty}^{+\infty} f(\mathbf{x})\ell(t,\mathbf{x})dm(\mathbf{x}).$$

for every regular function *f* , where (*t*, *x*) = *<sup>t</sup>* <sup>0</sup> 1*∂*O(*Xs*(*x*))*d*(*s*). Then we can consider stochastic differential equations in domains O ⊂ <sup>R</sup>*d*, of the form

$$dX\_t = \sigma(X\_t) \, dW(t) + u(t, X\_t)dt + n(X\_t)d\ell(t), \quad X\_0 = x\_{\prime t}$$

where *n* denotes the unitary vector normal to the boundary. These stochastic processes are reflected at the boundary and, in terms of partial differential equations (their generators), they correspond to considering Neumann boundary conditions.

Finally we consider the more recent notion of *forward–backward stochastic differential equation* (cf., for example [49]).

For given (smooth) coefficients *σ*, *u*, *v*, initial condition *X*<sup>0</sup> and final condition *h*, one looks for semimartingales *Xt*,*Yt*, *t* ∈ [0, *T*] which are solutions of the following system of stochastic differential equations (written here in integral form),

$$\begin{cases} \begin{array}{c} X\_t = X\_0 + \int\_0^t \sigma(X\_s) \, d\mathcal{W}\_\mathbf{s} + \int\_0^t \mu(\mathbf{s}, X\_\mathbf{s}) ds \\\ Y\_t = h(X\_T) - \int\_t^T Z\_\mathbf{s} \, d\mathcal{W}\_\mathbf{s} + \int\_t^T \nu(\mathbf{s}, X\_\mathbf{s}, Y\_\mathbf{s}, Z\_\mathbf{s}) ds \end{array} \end{cases}$$

with *E <sup>T</sup>* <sup>0</sup> [|*Xt*| <sup>2</sup> + |*Yt*| <sup>2</sup> + |*Zt*| <sup>2</sup>]*dt* < ∞. There are two remarkable features of such systems: one is that, in spite of a condition given at a final time (*h*), these solutions turn out, in fact, to be adapted to the past filtration. Another one is that, even if a priori we have three unknowns *X*,*Y* and *Z*, the last one will be in fact equal to *Zt* = ∇*u*(*Xt*). These kind of systems are natural generalizations of second order ordinary differential equations to the stochastic setting.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
