Effect of Transport Noise on Kelvin-Helmholtz instability

The effect of transport noise on a 2D fluid may depend on the space-scale of the noise. We investigate numerically the dissipation properties of very small-scale transport noise. As a test problem we consider the Kelvin-Helmholtz instability and we compare the inviscid case, the viscous one, both without noise, and the inviscid case perturbed by transport noise. We observe a partial similarity with the viscous case, namely a delay of the instability.


Introduction
Stochastic transport is a new fundamental perspective on fluid dynamics, see e.g. [18,23], [27] and [14]. A transport type noise in a fluid dynamic model may be seen, loosely speaking, as a simplified description of small-medium spacescales of motion. In the numerical simulations (see for instance [27]) we may observe the way it perturbes large scale motion; in general, this perturbation destabilizes large scales producing smaller eddies.
In this note we want to explore a special way transport noise may affect large scales, somewhat opposite to the one mentioned above. The key difference is the assumption that it is very-small-space-scale. The noise used below is made of very small, low intensity, vortex structures. In such a case it may happen that the transport noise acts as a dissipation, an additional viscosity.
It corresponds to Joseph Boussinesq intuition [5] that "turbulent small scales may be be dissipative on the mean flow". The physical intuition, beyond the specific mathematical derivation, is that fluid particles move so erratically to produce effects similar to the molecular motion. Although a proof is missing and the empirical validity is only moderate [24], this idea is certainly very useful for numerical simulations of fluid dynamics model, being the basis of the LES method [4], and, in case of vorticity equations, the vortex blob method [7].
Theoretically, this kind of turbulent-small-scale transport noise has been investigated in rigorous works. This line of research was initiated in [16] and developed in several works, see e.g [10,12] and many others (see [14]).
However, its stabilizing power has not been tested numerically yet. Here we observe its action on one of the strongest and most common instabilities: the Kelvin-Helmoltz one. We compare the inviscid case, the viscous one, both without transport noise, and the inviscid case perturbed by transport noise. The results are described in Section 4. Since we approximate the fluid dynamic equations by vortex methods, which fits quite well and in a unified way with the three cases analyzed here (inviscid, viscous and stochastic transport), we describe some preliminaries on this topic in Sections 2 and 3.

Model formulation
We consider the two-dimensional flow of an incompressible fluid on a 2D domain, either the full plane, T 2 or S 2 . As usual, the equations for the motion of the fluid are the conservation of mass and linear momentum, see e.g. [2], expressed through the Navier-Stokes equations in the null-divergence formulation. The main interest for our numerical simulations is the equation for the evolution of vorticity ω(t, x), which can be derived from the Navier-Stokes equations: Here ν is the kinematic viscosity, and u is the velocity of the fluid solving the NS-equation. In the 2D case, we can express ω := ∇ × u as the curl of the velocity field, a vector normal to the flow plane. We denote x = (x 1 , x 2 ) to be any point in the domain.

Point vortex method for inviscid flows
In order to state the problem we investigate numerically, we first recall several well-known facts; for a general introduction, see e.g. [21,22,26]. For the sake of simplicity, we focus on an inviscid fluid; the vorticity equation (1) reduces to, From the null-divergence hypothesis, we define the stream function associated to the fluid ψ(t, x), see [2], such that the velocity u is given by u = −∇ ⊥ ψ. We obtain the stream function by solving the Poisson equation To compute the solution, we need to express the Green's function and the convolution with the variations of constants methods, which for flows in the full plane R 2 reads as where G is the Green's function, the fundamental solution of the Laplace equation. The explicit form of G in the full plane R 2 is Using (4) we obtain the velocity field where K is given by The equations (6) and (7) for the velocity are known as the Biot-Savart law, and K is the Biot-Savart Kernel. Note that we must correct this velocity for flows over bounded or periodic boundaries domains to satisfy the boundary condition. Therefore, we have to change the Green's function according to the Poisson equation. Consider now a fluid particle X t , moving in the velocity field; from (2), the path of the particle is Therefore, considering the fluid as distinct "fluid particles" of constant vorticity, these abstract objects' motion determines the scalar field's evolution. This is the premise of the point vortex method to compute (2), (see [21,22] and [26] for details).
To this end, consider N point vortices, idealizing a 2D inviscid fluid and occupying positions X 1 t , ..., X N t , with intensities (circulations) Γ 1 , ..., Γ N respectively. They move accordingly to the following set of ordinary differential equations derived from the previous computations (8): where the vector-valued kernel K (x, y) is the Biot-Savart kernel, equal to 1 2π (x−y) ⊥ |x−y| 2 in full space, suitably modified on a torus or in a bounded domain. One can prove that the empirical measure is a weak solution of 2D Euler equations in vorticity form (suitably interpreted for distributional fields as in [25]): with Here ω is the (scalar) vorticity, u is the (vector) velocity.
Given a bounded probability density ω 0 , taken a sequence of i.i.d. r.v. X i 0 with density ω 0 , taken Γ i = 1 N above, considered now ω 0 as a random variable with values in distributions, the random empirical measure converges weakly, in probability, to the unique solution of the above Euler equations (9) (when the initial condition ω 0 is measurable and bounded, the Euler equations have global existence and uniqueness of bounded measurable solutions). Similar results hold also when the approximation of the initial condition is deterministic, and in the case when ω 0 is more general than a probability measure, with suitable modifications of the scheme. Notice that the velocity field associated to the distributional vorticity ω N (t, ·) is This is a well defined vector field, of class L p loc for every p < 2 but not for p = 2.

Point vortex method for viscous flows
To investigate viscous flows, we modifying the previous scheme by adding independent 2D Brownian motions W 1 t , ..., W N t to the equations of point vortices Then, the empirical measure ω N (t, ·) converges weakly, in probability, to the unique solution of the 2D Navier-Stokes equations in vorticity form 3 Point vortex method with environmental noise As mentioned above in the Introduction, several works indicated an interest in the following stochastic modification of the Euler equations where σ k = σ k (x) are given vector fields, that we assume divergence-free, B k t k∈K are independent 1D Brownian motions and the stochastic operation • stands for the Stratonovich integral. Due to this, formally, vorticity is conserved (it is transported randomly by the field udt + k∈K σ k dB k t ). This model bears similarities with the viscous flows, and the objective of this paper is to show differences and similarities of the elliptic operator obtained from such a transport-advection noise.

Transport noise and deterministic scaling limit
The point vortex dynamics associated to the model (12) is then given by the following expression: Notice that this is a model of common noise (also called environmental noise): the BM's B k t are the same for all particles, opposite to the model (11) where each particle X i t was affected by an independent BM W i t . See [11] for an example of theoretical results on this model. For models similar to this one, it has been proved (see e.g. [8]) that the empirical measure converges to the solution of the SPDE (12). At the same time, following [16] and subsequent works, if the noise is parametrized in such a way to become more and more small scale, the SPDE (12) converges to the deterministic equation with additional viscosity Inspired by [13], we consider a sort of mixed scaling limit: we take the point vortex dynamics with common noise (13), which for given fields σ k would converge to the SPDE (12), and we choose more and more small scale coefficients σ k in order to be close to the deterministic equation (14). The present work is numerical, but the theoretical scaling limit behind it would be that the point vortex model (13) converges to the Navier-Stokes equation (14). Recalling the result mentioned in the previous section, namely that point vortices perturbed by independent BM's (10) also converge to the Navier-Stokes equation (14), we see that two different models of noise, (13) and (14), lead to the same limit equation. Our aim is to explore the validity of this fact from a numerical viewpoint and in the particular case when also the coefficients σ k are point vortices.

A digression on the theoretical selection of the noise
In this section, in the same spirit as [13,14,15], we explore some property of the environmental noise that we use in a simplified way in our numerical simulations. Following the works on modeling of passive scalars [19], when considering the scaling limit of [13] to ω(x) solution of the viscous Euler equation [14], we consider a model of noise in the fluid which is delta-correlated in time, namely a white noise with a precise space dependence.
where (σ k (x)) k is a family of smooth divergence free vector fields on the 2D domain of the equation, and B k t are independent one-dimensional Brownian motions; K is, usually, a finite index set, but with suitable assumption we could consider also the case of countable family of smooth fields. In this case, the term W (t, x) · ∇ω (x) obtained in the convergence result of the point vortex empirical measure, must be interpreted as a Stratonovich integral This is given by an Itô-Stratonovich corrector plus an Itô integral; precisely, is given by: where M (t, x) is a (local) martingale. Follows that the Itô-Stratonovich corrector takes the form of an elliptic operator: As an example, we take the noise [19], which is relevant to our numerical investigation in the choice of the divergence-free field in the point vortex model. For simplicity, assume the domain to be R 2 , but modifications on T 2 , S 2 are possible. Its covariance function is space-homogeneous, i.e. C (x, y) = C (x − y), with the form The famous Kolmogorov 41 case follows if we take ζ = 4/3. Taking k 1 = +∞, then C (0) = Kσ 2 where the constant K is given by We consider small-scale turbulent velocity fields depending on a scaling parameter and taking the scaling limit in 12, as in [13,16]. In the case of [19] we have and simultaneously, we may have that the Itô term goes to zero, hence recovering 14. Since this procedure can be imagined as a sequence of scaling limits with K, N being large, we can't expect a precise convergence to the Laplacian operator at the level of a numerical study, but a similar effect is expected on the fluid, namely: a diffusive behaviour and a delayed formation of classical pattern and large scale structure in the fluid vorticity.
4 Numerical results

Setting: Kelvin-Helmholtz instability
In this section, we investigate classical results on the shear flow model in the setting of point vortices, analyzing the Kelvin-Helmholtz instability and the possibility of delaying the structure formation. In this way, we can both test the goodness of our point vortex models and, at the same time, set a benchmark for which we will show delayed instability. In order to test the point vortex model in the classical cases (8) and (10), we choose the particular fluid configuration of a shear flow because of its fundamental property: developing instability without viscosity and delaying it when viscosity is present. We work on a strip S 2 equal to the set [−1, 1] × R with coordinates x = (x 1 , x 2 ) and identified boundaries at x 1 = ±1; all fields are periodic in the x 1 -direction. We take an initial velocity u 0 of the form and corresponding vorticity ω 0 = ∂ x2 u 0 1 (x 2 ). We choose, in particular, the function where we fix δ = 0.02 in our numerical simulations.
To compute the vorticity measure of our point vortices, we use vortex blobs, obtained by spreading the circulation of a point vortex over a chosen small area, the vortex core (see e.g. [26]). In this formulation, the vorticity field is approximated by where the mollifier φ ε describes the vorticity distribution in the vortex core, the subscript ε represents the characteristic size of the vortex core. Following standard numerical techniques (see e.g. [3]), the core size ε of the vortices has to be much larger than the average spacing d between the vortices; the core size is usually taken to be ε = d q , with q << 1.
In (17), the vorticity distribution at any time depends on the point vortices X i t through the vortex blobs. In our numerical simulations, we take N ∼ 10 4 point vortices; following a mean-field approach, the initial circulation for every given point vortex is derived from u 0 1 and is equal to Γ i 0 = 1 2δN . We solved the point vortex model (8) using Heun's method for a second-order time discrete approximation; the time step for our simulations was selected to be ∆t ∼ 10 −3 to ensure a trade-off between the stability of our method and the generation of vortex-like structures in the shear flow model.
In the usual way, we also recall that in our numerical framework, the kernel K in (13) corresponds to the Biot-Savart kernel. We have that where G is the Green function on S 2 . In the whole plane we have the simple expression G R 2 = 1 2π log |x|; while for our domain we know that and s(x) is a smooth function on S 2 . Thus, K is divergence-free, smooth away from the origin, and symmetrical; moreover, it holds the following behaviour: |K(x)| ∼ 1/|x|, as |x| → 0 , which we extensively use to approximate our kernel with K R 2 , throughout the numerical simulations. Without ambiguity, from here on, we consider the horizontal and vertical axes as our reference frame, naming them the x-axis and y-axis, as usual.

The role of intrinsic instability
We know that, at least formally, the vector field u = u 0 is a solution of Euler equation (14) with ν = 0. This system is unstable: small perturbations rapidly develop vortex blobs. As such, we consider the system of point vortices (X i t ) i with initial vorticity expressed as where the circulation for each point is equal to 1 2δN and the initial positions of the vortices X i 0 ∀i = 1, ..., N are uniformly distributed on the strip [−1, 1] × [−δ, δ]. The randomly generated initial condition represents small perturbations in the system and is responsible for the different pattern formation.
The measure ω N t := 1 N i δ X i t converges, in distribution, to the scalar vorticity field solution of the Euler equation (14); analogously with the continuous The exact solution of the Navier-Stokes equation (14), with ν > 0 and initial condition u 0 , is given by where u 1 (t, x 2 ) solves the heat equation, Due to the spreading of the profile u 0 1 , the solution becomes more stable; namely, the development of vortex blobs is delayed. In our numerical simulations, we reproduce this phenomenon by perturbing the system 8 through independent Brownian motions B i t , i = 1, ..., N , with variance linked to the viscosity parameter: V ar(B i t ) ∼ √ ν. This system converges to the exact solution when N → ∞; however, in our numerical study we are dealing with a finite system. For this reason, the profile of the strip remains quite stable for short times, with just a spread along the y-axes. We report in figures 2a and 2b a single configuration at two different timesteps, t = 50 and t = 100; we take √ ν = 0.095 to better focus on the stability restoration. Our results are in agreement with the theory: from a comparison with 1b-1c, we see that when ν > 0, the profile is much more stable and diffused than in the deterministic case, and blob-like structures appear only at large times.

Selection of divergence free field
Starting with the same initial condition (16), we consider N + M point vortices, each of which we associate with a position in S 2 in the following way: Here, the vortices Y i , i = 1, ..., M do not move, and when activated, they represent the feedback of small-scale turbulence acting on the fluid itself on large scales. The new simulated vortex dynamics for X i t as in (13), reads where the Brownian motions are all independent, B i t bi-dimensional and W j t uni-dimensional, and they are acting simultaneously on all the particles i = 1, ..., N . The environmental noise follows the Stratonovich integral prescription, automatically implemented in Heun's method [20]. We choose the diverge-free vector fields σ j as following the theoretical analysis performed in [14,13]. Here, the intensities a N,M j are linked to the scaling limit process, which produces a viscosity term on the large scales, and K, the Biot-Savart kernel, simulates the action of such small vortices. The idea behind such a selection is that we want to exploit the same features of the vortex model, with the formation of small-scale vortex structures generating feedback on the entire configuration. In the limit, the dynamics of such small structures, modulated through a Brownian motion, rebound on large scales, perturbing their motion with their dissipative properties and delaying the formation of the instability.

Positions and intensities of fixed vortices
In order to make contact with previous studies [13], we choose the positions of the fixed vortices Y j and their intensity a N,M j according to the convergence of the scaling limit [4.2.1]. More precisely, at each timestep, we generate Y j , j = 1, ..., M , uniformly distributed point vortices; their position on the y-axis is apriori selected in the interval [−δ F X , δ F X ]. In this setup, the vortices Y j are generated in a strip of variable height 2δ F X ; this strip contains the moving vortices X i t , and it is taken to be of the same height of our boundary fluid layers, or one order of magnitude greater. This choice emphasizes that our proposed "small-scale" structures should act on all points of the fluid in all directions: the average contribution of the Y j on the X i t along every direction should mimic a Brownian motion. We explored different setups of positions and intensity; we selected meaningful realizations, as reported in table 1.  We choose the intensity of the "small-scale" perturbations following heuristic considerations. We consider the mean inter-particle distance between two fixed point vortices, < r >= √ m, where m = A/M is the particle density and A is the total area occupied by the M vortices.
Then, we focus on a single moving vortex, X i t ; we compute the magnitude of its velocity when X i t is at a distance d =< r > /2 from the nearest fixed vortex, i.e., its position is halfway between two fixed vortices. As a result, we obtain the following estimate for the velocity of X i t : where we suppose that the coefficients depend on the configuration (X i t ) i , and are equal for each j = 1, ..., M . What is left is to estimate the number of fixed vortices such that the interaction with X i t is not negligible: let us call this number K, giving us ∼ Ka N 4πd . Thus, using the theory from the scaling limit of environmental transport noise (see e.g. [14,13,16]) and the construction of section 3.2 for point vortices with transport noise, we assume that This leads us to the estimate for the intensity of the fixed vortices It remains to estimate K, the number of the nearest fixed vortices: consider a ball centered in X i t with radius d, so that the area is A near = πd 2 . We recall that m is the density of the fixed vortices, then the nearest vortices are: Taking into account only the nearest vortices, we are underestimating the actual contribution of all the vortices. In particular, we should compute such contribution by considering a radius dependent on the range of the image of the Biot-Savart kernel. For this reason, we empirically selected a wider radius αd, with α ∼ 3. Concluding we get our estimates for the intensity

Effect of small scale common noise
As recalled in the previous paragraph's heuristics, the procedure of the scaling limit is preserved when both N and M are large, and the intensity a N,M j is small. For this reason, we do not search for the same exact solution of the Navier-Stokes equation (14), ν > 0, with initial condition u 0 , as per the case of the independent noise. However, since the regime tends, in the limit, to the same solution, we expect a diffusive effect on the strip of the point vortex. More precisely, we expect to see a delay in the formation of macroscopic structures and a more dispersed displacement of such small vortex blobs that, on average, should maintain the strip configuration for a longer time.
In the first of our simulations, we generate at each time step M ∼ 2 · 10 5 fixed vortices, with intensity a N,M j ∼ 5 · 10 −4 which follows from the heuristics. The fixed vortices are uniformly distributed in a strip [−1, 1] × [−0.1, 0.1], containing the initial point vortices configuration. This particular setup captures the feedback effect of small scales on large vorticity structures, as the contribution of all the low-intensity perturbations on the dynamics of the point vortices averages in every direction. In the proposed numerical simulation, we show that the transport noise model reproduces the desired instability delay, even if it is slightly less effective than in the independent noise case; we illustrate a snapshot of a configuration for time t = 50 and t = 100 in figures [3a,3b].
From a comparison with figures [1b,2a], we see that the initial strip configuration is preserved for a longer time than in the deterministic case and rotation of the fluid is milder, but the profile is less stable than in the viscous regime. If we focus on the deterministic case, we see blob-like structures formation already at t = 50; in contrast, in the transport noise regime, those structures are less visible and appear more prominently only at the end of our simulation (t = 100). This delay of the instability is evident in the realizations in figures [3b], compared with [1c,2b]: we notice a more diffused and homogeneous profile and a delayed formation of rotational structures due to the noise spreading the particles along the y-axis. A difference with the viscous case is that the compression in the x-axis is stronger than in the case of the independent noise, resulting in a more prominent stretch, which could resemble more the deterministic formations, placing the transport noise as a midpoint between the two regimes.
In the second of our simulations, the strip of fixed points Y j is generated in the same region as the point vortices at each timestep; we selected M ∼ 1, 32 · 10 5 , the fixed vortices are uniformly distributed in [−1, 1] × [−δ − , δ + ], with = 0.03, and their intensity is derived from the heuristics a N,M j ∼ 5 · 10 −4 . The results are shown in figure [4]: diffusion on the y-axis is present for short times, and preservation of the strip profile is guaranteed. However, the drawback of such a configuration is that analysis can be performed only for a short time: boundary effects of the fixed vortex strip can deteriorate the configuration, making the results unrealistic. In future works, we expect to overcome this obstacle by proposing a new method, now in the study, to generate small vortices only in regions activated by the shear flow's movement. We performed a final simulation, in which we take the density of fixed vortices to be smaller than the density of the point vortices. In particular, the strip

Diagnostics
In this section, we perform stochastic analysis on the three configurations proposed in this study to highlight the differences and the reconstructed stability, or the emergence of new structures, in the Kelvin-Helmholtz instability problem. As a first step, we compute the vorticity ω obtained through the same mollification as in Majda [3], through the vortex blob method applied to each of the point vortices X i t . We report our results for the vorticity computed in the deterministic case at t = 100 in figure [6a]. We see that the vorticity measure concentration near the fluid's boundary layer is located in the newly developed macroscopic structures. Moreover, a displacement from the initial configuration where the laminar fluid started its evolution is present.
In contrast, the vortex blob solution with viscosity ν > 0 retains its structure for longer times than in the inviscid case. The instability delay is graphically evident both from the configuration reported in figure [2a] and the vorticity intensity reported in figure [6b]: the vorticity measure concentration at time t = 100 is similar to the one of the initial strip but with a more diffused profile on the horizontal line.
Finally, in figure [6c], we show the vorticity in the environmental noise regime at t = 100. In contrast with the inviscid case, we see no development of macro-scopic structures in the profile. However, even though a more diffused profile, with lower density overall, is present, the stability is lost at larger times. This instability at larger times suggests that different behavior, dependent on the density of fixed vortex Y j and selection of transport noise fields σ j , could arise in applying this kind of small-scale approximation. For this reason, we focus our analysis on small-time behaviour. Concerning the formation of large rotating structures, we see that the particles spread in the horizontal direction when a forcing term, either an independent or transport noise, acts on the fluid, in contrast to the solution of the Euler equation with ν = 0. In particular, we focus on the empirical density obtained from the x-axis in the three configurations at t = 100, [7]. The deterministic case [7a] shows a complete formation of separate blobs with peaks in the exact locations of the macroscopic structures, as in 1c. On the contrary, in the viscous 7b and transport noise case 7c, the distribution of the vortices is more uniform, and it delays the instability of the fluid layers.
Following theoretical results, we know that with our initial condition ω 0 , the velocity u t solves 4.1.2 when ν > 0. This is crucial to understanding our system's short-time behaviour and the preservation of the initial configuration. This result states that the empirical density obtained from the y-position, when viscosity is present, maintain a Gaussian profile through time. In particular, in the case of ν = 0, the viscosity follows a classic Euler equation. As such, the yposition profile is far from a Gaussian: it behaves like a multi-modal distribution, concentrated in the proximity of the center of the large structures. This result is supported both by the profile of the particle system and by figures [8a,9a]; the qq-plot shows a distinct behaviour for small quantiles. The Kolmogorv-Smirnov test estimates a D-statistic of 0.031, with a p-value less than 10 −9 confirming the rejection of the Gaussianity hypothesis. Vice versa, when viscosity is present, i.e. ν > 0, as in the case of independent Brownian motions, the noise's diffusive behaviour allows the profile's restoration in the y-direction: the strip configuration is preserved for a longer time. From the profile of the point vortex system and figures [8b,9b], we see that the empirical density approximates well the one of a Gaussian kernel. Moreover, the qq-plot suggests a perfect match with a normal distribution, suggesting the preservation of the strip through time, trading it with more spread on the y-axis. Finally, performing a Kolmogorv-Smirnov test, we see, in fact, a D-statistic of 0.005 and a p-value greater than 0.9 suggesting to accept the normality hypothesis. This behaviour is preserved throughout the simulation, degrading only at longer times when a few large formations start to rise, as in figure [2b].
In the transport noise case, we know that we recover the same viscous Euler equation solved in the independent Brownian motion case by applying a particular scaling limit procedure. To this end, we expected that the more stable profile shown in figure [3a] presents the same diffusion on the density of the y-position as the case of the independent noise. This is the case as reported in figures [8c,9c]: the behaviour at a short times, t = 50, in which the profile still approx- imates a Gaussian kernel. We perform a qq-plot and a Kolmogorov-Smirnov test on the y-position; we obtain a D-statistic of 0.011, and a p-value of 0.312; those results suggest the profile stability and the validity of our hypothesis. However, later on (t = 100), even though the quantity obtained from the KS test is still preserved as in the viscous case, with a D-statistics of 0.008 and a p-value of 0.48, the profile degrades as shown in [3b]. Those results show that the strip configuration is non-preserved for longer times due to transport noise stretching acting on the point vortices. This can be seen in the tail of the distribution in figure [9c], compared to the viscous case in [9b], which shows at t = 50 already a different behaviour. This suggests that the environmental noise's ef-fect is not only related to the diffusivity of the strip, but is also responsible for the stretching and formation of different structures.

Concluding remarks
In the present work, we have produced numerical simulations of 2D incompressible fluids, also perturbed by transport noise, using the point vortex method. We focused on the special case of shear flow formation that produced a Kelvin-Helmholtz instability, in order to test the dissipativity properties of small-spacescale transport noise. We confronted the intrinsic instability generated in the deterministic case with the possible recovery of the stability through injected noise in the system in the form of transport noise. We showed that, for short times, with a degree less intense than the viscous case, we can maintain the stability of the strip at the expense of a more small-scale irregularity and diffusion of the profile.