Partitioning Waves and Eddies in Stably Stratified Turbulence

: We consider the separation of motion related to internal gravity waves and eddy dynamics in stably stratiﬁed ﬂows obtained by direct numerical simulations. The waves’ dispersion relation links their angle of propagation to the vertical θ , to their frequency ω , so that two methods are used for characterizing wave-related motion: (a) the concentration of kinetic energy density in the ( θ , ω ) map along the dispersion relation curve; and (b) a direct computation of two-point two-time velocity correlations via a four-dimensional Fourier transform, permitting to extract wave-related space-time coherence. The second method is more computationally demanding than the ﬁrst. In canonical ﬂows with linear kinematics produced by space-localized harmonic forcing, we observe the pattern of the waves in physical space and the corresponding concentration curve of energy in the ( θ , ω ) plane. We show from a simple laminar ﬂow that the curve characterizing the presence of waves is distorted differently in the presence of a background convective mean velocity, either uniform or varying in space, and also when the forcing source is moving. By generalizing the observation from laminar ﬂow to turbulent ﬂow, this permits categorizing the energy concentration pattern of the waves in complex ﬂows, thus enabling the identiﬁcation of wave-related motion in a general turbulent ﬂow with stable stratiﬁcation. The advanced method (b) is ﬁnally used to compute the wave-eddy partition in the velocity–buoyancy ﬁelds of direct numerical simulations of stably stratiﬁed turbulence. In particular, we use this splitting in statistics as varied as horizontal and vertical kinetic energy, as well as two-point velocity and buoyancy spectra.


Introduction
Most geophysical flows on Earth exhibit density variations and are submitted to the effect of rotation, so that their dynamics is the result of several intermingled phenomena of different nature: turbulent transport, waves, thermal diffusion, and convection, to cite but a few. Focusing on density variations, stratification is present, e.g., in convective layers wherever the situation is unstable-light fluid toped by heavy fluid-as in the troposphere, or in a more stable setting in the stratosphere. These effects are mostly due to potential temperature gradients in the atmosphere, but also appear at different scales in the ocean, due to salt-and temperature-related density variations [1]. In a stably stratified flow, with negative density gradient aligned with the gravity axis, one observes the presence of internal gravity waves (IGW) at different scales and frequencies. They therefore contribute to the flow dynamics, in terms of transport properties, and of interscale energy transfer. Moreover, in investigations show that this technique overestimates the wave energy as it does not take into account the temporal dependence [23]. In addition, it does not directly apply to other flow types such as rotating or rotating-stratified turbulence.
We therefore seek additional criteria for evaluating the presence of waves in a flow, and we wish to propose further insight into the analysis of wave-related dynamics in turbulence. As befits the definition of waves via their space-time coherence, the dispersion relation that defines inertial gravity waves in a turbulent flow links a spatial component-the angle θ-to a temporal component-the angular frequency ω-and can be applied to any wave that satisfy a dispersion relation.
Before considering in a further study the general case of rotating and stratified turbulence, we propose in the present study to consider theoretically and numerically how waves in stratified turbulence are created and how they are convected by simple large-scale flows whose complexity is progressively increased. In the following, we first apply a straightforward spatiotemporal statistical analysis to different canonical flows whose dynamics is assumed to be linear, and we observe that the sweeping effect appears to be the main source of wave frequency modification. Then, a new technique based on spatiotemporal decomposition is implemented that requires a four-dimensional (4D) Fourier transform (3D in space and 1D in time) and permits to filter the waves by using the dispersion relation modified by the sweeping effect due to VSHF. With this decomposition, statistics for the separate waves and eddy contents of the flow are obtained, e.g., their energy spectra, as well as maps of the respective flow structures by visualization in physical space.
In the following, we present the equations and statistical flow characterization in Section 2, the results of the sweeping and Doppler canonical effects in Section 3 and our spatiotemporal decomposition with detailed analysis and comparison of wave-vortex flow decomposition in Section 4, before concluding with general remarks and perspectives in Section 5.

Governing Equations and Dispersion Relation
We consider the Navier-Stokes equations in the Boussinesq approximation for an incompressible fluid: where u = (u x , u y , u z ) is the velocity vector, ω = ∇ × u is the vorticity, b is the buoyancy field, p is the pressure field, F u and F b are the external forcing, ν is the kinematic viscosity, and χ is the thermal diffusivity. n is a vertical unit vector for the axis that bears gravity (along the vertical direction z). The buoyancy b accounts for the fluctuation of density around the mean gradient N 2 where N is the Brunt-Väisälä frequency. For simplification, the Prandtl number is chosen equal to one: Pr = ν/χ = 1. The fields in these equations are assumed to be non-dimensionalized by reference length and time scales. A plane wave solution exp(i(k · x + ωt)) with wavevector k and frequency ω is obtained for the linearized Equation (1) without viscosity. Viscous dissipation can easily be re-introduced, although the inviscid solution is retained here for the sake of simplicity. From this inviscid plane-wave solution, the dispersion relation of internal gravity waves [5] can be recovered as ω r = ±N cos θ (2) where θ is the angle of the wavevector k with the horizontal plane (x, y) (see Figure 1a). This dispersion relation is very particular because it involves the stratification intensity N and the angle of propagation θ of the wave. Furthermore, unlike sound waves or surface waves, the dispersion relation depends only on the orientation of the wavevector k = (k x , k y , k z ) and not on its amplitude k = |k|. We denote k h = (k x , k y ) the projection of k on the horizontal plane (see Figure 1a), and k h = |k h |. According to Pedlosky [24], the crest of constant phase φ = k · x + ωt = constant propagates perpendicularly to the wavevector k since ∇ x φ = k with a phase velocity v Φ = ωk/k 2 . It propagates with a phase velocity v Φ parallel to k if ω < 0 and anti-parallel to k if ω > 0. In Figure 1b, we plot the four possibilities (the four quadrants Q 1 , Q 2 , Q 3 , Q 4 ) of the propagation of the crest of phase in the plane (θ, ω) according to the signs of θ and ω. Note that the envelope of the wave, (i.e., the energy of the waves) propagates along the direction of the group velocity v g = ∇ k ω which is perpendicular to the wavevector k for IGW [24]. To generalize our analysis, we introduce k = (−k x , −k y , k z ) as the horizontal opposite of k and sharing the same angle θ.

Numerical Method
The complete nonlinear Equation (1) are non-dimensionalized and solved using a classical pseudo-spectral method [25] in a periodic domain of length 2π in each direction. The numerical domain is periodic with a grid spacing uniform and a number of grid points noted n g in the three directions (n 3 g points overall). The 3D spatial Fourier transform is done with the P3DFFT library [26]. The viscous and diffusive terms are solved via an implicit treatment implemented as a change of variable involving an exponential integrating factor. The transformed variables for the velocity vector and buoyancy fields are in Fourier space: thus reducing a possible numerical approximation error for diffusion. The dealiasing of the quadratic terms is based on the phase shifting method [27] with a truncation at k max = √ 2(2/3)k N with the Nyquist wavenumber k N = n g /2. A third-order Adams-Bashforth time-scheme is used to perform numerical time integration. The code is parallelized using MPI. Simulations are done either in the laminar case or in the turbulent case. In the laminar case, a specific forcing F u and F b is chosen which is explained case-by-case in the next sections. In the turbulent case, a specific spectral statistical forcing (explained in [28]) is used, with the advantage of limiting the increase of VSHF. It consists in forcing velocity only (F b = 0) with a forcing term F u chosen in Fourier space to ensure a constant injected power P =< F u · u >, where < * > is a statistical volume average. The force F u is applied in a cylindrical shell of radius k v with finite vertical extent between k v,min and k v,max . The mesh size is chosen to ensure accurate resolution of the smallest scales through the condition k max η 1, where k max is the maximum resolved wavenumber, η = (ν 3 /ε) 1/4 is the Kolmogorov length scale and ε is the kinetic energy dissipation. The space-time evolution of the velocity and density fields are then post-processed to obtain various statistics, as in Section 2.3.

Space-Time Analysis of the Flow
Space-time statistics are computed to link the direction of propagation θ-a spatial variable-to the frequency ω-a time variable. This is done by calculating the time-dependent angular-dependent spectral density of kinetic energy. To do this, the wavevector space is first decomposed into different elementary cones with an angle θ with respect to the horizontal plane and a length k (see Figure 2). Due to the θ angle discretization, the elementary cone is a thin volume where all wavevectors k have an angle θ(k) ∈ I θ = [θ − ∆θ, θ − ∆θ]. To compute the energy into this elementary cone, we decompose the computation into four steps: 1. Each elementary cone, which contains wavevector k with an angle θ(k) ∈ I θ , is divided into sub-volumes by selecting the scale k = |k|. Then, all the velocity vectors are summed onto the sub-volume in a new variable U(θ, t, k) = ∑ θ(k)∈I θ ,|k|=k u(k, t).

The Fourier transform in time of that new variable is computed asŨ
3. The spectral density of kinetic energy is recovered in a (θ, ω) plane for a range of scales between k 1 and k 2 that is a set of sub-volumes: 4. The total spectral density of kinetic energy into an elementary cone, regardless of scale k, can be computed as a summation over all sub-volumes: When the numerical simulation contains only waves, one should observe a peak of energy density that follows the dispersion relation curve defined by Equation (2) in the (θ, ω) plane (as illustrated in Figure 3b). We detail this in Section 2.4. Figure 2. Decomposition of the Fourier space in elementary "cones" containing wavevectors k with given wavenumber amplitude k, and θ within discretized intervals between 0 (k z = 0) and π/2 (k h = 0).

Saint Andrews cross Wave Propagation Benchmark
We test here the statistical characterization of waves by a simple benchmark consisting of a single wave propagation in a Saint Andrew cross pattern, as in experiments [5]. The flow forcing produces only waves and is examined from its initial condition at rest. Several forcing possibilities exist to create either a single isolated wave or several superimposed ones. The forcing term is introduced in the linearized Navier-Stokes equations in the Boussinesq approximation: the non-linear terms ω × u and u · ∇b are removed from Equation (1), likewise in the numerical simulation. For this benchmark, the forcing F b is only imposed in the buoyancy equation, so that F u = 0. We choose a point forcing localized in physical space, of the form: with ω f = 0.3. This forcing implemented only in the buoyancy equation of system Equation (1) ensures that the incompressibility condition is maintained for the velocity field. The function δ x is the Dirac function. In the Fourier space, it means that all wavenumbers are forced even though there is more energy for high wavenumbers. We use this forcing in a simulation without non-linear terms. The resolution is n g = 128, we assume negligible viscous diffusion so that ν = 10 −7 , and the Brunt-Väisälä frequency is N = 1. As expected, we observe in Figure 3 the propagation of waves according to an angle θ set by ω f such that ω f = N sin θ. Figure 3a shows the distribution of the vertical velocity components, and the same pattern is observed in the buoyancy field (not shown). The spatiotemporal analysis is then applied using 1000 fields separated by a time step ∆t = 0.5. The total time span is therefore 79.6T N where the T N = 2π/N is the Brunt-Väisälä period. Upon computing the space-time statistics of energy density E(θ, ω) according to Equation (5), the energy distribution in the (θ, ω) plane (Figure 3b) appears to concentrate in two kinds of regions: (a) along two horizontal lines of energy concentration, at ω f = ±0.3 (these lines are the traces in the numerical simulation of the pointwise spatial forcing-almost a Dirac-that oscillates at frequency ω f , which is Fourier transformed); and (b) along the dispersion relation curves (in absence of non-linearity, no energy redistribution can occur away from the immediate input due to the forcing). Both kinds of energy concentration curves correspond to the response of the linear system to the forcing, with one component at its forcing frequency ω f and one at the natural frequency of the system, which is here the frequency of the waves ω r . The maximum of E(θ, ω) is therefore observed at the crossing points of the forcing and the dispersion relation curves, that is at θ 72.5 • . In physical space, this peak energy results in the observation in of waves propagating at θ 72.5 • , which is indeed what we measure in Figure 3.
This numerically observed solution of the response of the linearized system of equations to harmonic forcing is computed analytically in Appendix B from the inviscid Navier-Stokes-Boussinesq equations with a zero initial condition. The analytical solution for the buoyancy field is: From this solution, two regions of concentration of energy can indeed be found: one along the line defined by ω = ±ω f , which is the frequency of the forcing, and one along the curve ω = ±ω r , which is the frequency of the waves along the dispersion relation. This solution diverges for ω f = ±ω r because it is an inviscid solution, but it shows that the maximum of energy is at the intersection of these frequencies.
Note that several vertical lines at constant angle θ are visible in Figure 3b in a log-scale. Each line based on constant θ corresponds to an independent set of Fourier transforms in time. We think that there are two main numerical reasons for these vertical lines. Firstly, numerically, the Dirac function in Equation (7) is replaced by a sinc function, with visible effect in log scale. Secondly, we do not use any windowing technique because it introduces an artificial coefficient that changes the energetic ratio between the dispersion relation and the input signal. This would make difficult to filter the waves from the eddies as done in Section 4. To remain consistent with the rest of the article, we do not apply a windowing technique (such as Hann) and we analyze the raw data. Nevertheless, in Appendix A, we apply this technique on the same data as in Figure 3b: this technique increases the concentration of energy around the dispersion relation. Moreover, this vertical line tends to disappear when the four-dimensional filtering is used (see Section 4.2) or when statistical stationarity is reached (see Section 4.6).

Sweeping and Doppler Effects on Internal Gravity Waves Propagation
We show above that the dispersion relation can be recovered accurately from a space-time statistical analysis of a simple flow where a single wave propagates, as shown in Figure 3. However, in more complex flows, waves can be convected by eddies, by a mean shear or by any other kind of large-scale flow. To understand how the frequency of the waves can be modified, and how this translates into a new pattern of energy E(θ, ω) concentration in (θ, ω) plane, we consider a mean flow homogeneous-i.e., periodic-in one spatial direction, and we compute its effect on internal gravity waves propagation.
The reader can find an overview of our simulation cases in Table 1 that gathers the main parameters for the laminar cases and the turbulent case. Results for the turbulent case are presented in Section 4.3. Table 1. List of numerical parameters for the different numerical simulations presented in this paper.

Sweeping by a Uniform Convective Flow
To begin with, we consider the case of a uniform convective flow with mean velocity vector either in the vertical direction-i.e., along the axis of propagation symmetry-or the horizontal direction-perpendicular to the axis of propagation symmetry. The sweeping effect of the flow is then characterized for both cases from both theoretical and numerical predictions.

Theoretical Modified Dispersion Relation
The advecting mean flow c = (c x , 0, c z ) is introduced in the Navier-Stokes-Boussinesq equations as follows: Since the mean flow is uniform, the convection term can be introduced in the numerical simulation via a change of variable (similar to that for diffusion in Equation (3)) where the new velocity vector and buoyancy fields are: With this new set of linearized equations, the inviscid dispersion relation becomes ω = ω + c · k, which shows that the wave deviation is larger at large wavenumber k. This new dispersion relation can easily be found by using a change of variable in the time and the spatial Fourier domains. According to Figure 1b, a kinematic analysis permits to infer the variation of the frequency in function of the advective flow c: when c is in the same direction as the phase velocity (see Figure 1b), the angular frequency increases, whereas, when c is in the opposite direction of the phase velocity, the angular frequency decreases. Such a modification of the dispersion relation was already observed in a rotating flow experiment, although indirectly: measurements of a Doppler-like frequency shift were explained by taking into account sweeping by the two-dimensional geostrophic mode observed in the turbulent fluid tank [29].

Sweeping by Vertical Velocity
We start by considering a vertical advecting velocity field c = (0, 0, c z ), i.e., along the natural axis of symmetry of the system, borne by the gravity vector. The Saint Andrew's cross pattern obtained without mean flow as in Section 2.4 is therefore convected towards positive z direction. As a result, the dispersion relation is modified by the term c · k = c z k z = c z k sin θ: The result is shown in the two panels of Figure 4 for c z = 10 −2 , with the same physical and numerical parameters as in Section 2.4. Figure 4a shows that the wave-packet is translated, thus producing a 'wake-like' pattern, much as that of a uniformly propelled boat for surface waves. Figure 4b shows the distribution of energy density E(θ, ω). This figure shows that the dispersion relation is modified according to Equation (10), and that maximum modification is observed when considering the largest wavenumber range 56 < k ≤ k max = 60. This domain of wavenumber is chosen in order to observe the maximum effect of the sweeping effect. Moreover, the energy density E(θ, ω, 56 < k ≤ 60) gets very large at the intersection between the forcing and the modified dispersion relation with an asymmetrical angle value of θ ±0.8 rad and θ ±1.3 rad. Note that Figure 4b is in agreement with Figure 1b: in quadrant Q 1 , c is in the opposite direction compared to the vertical component of phase velocity v φ , and the frequency ω decreases compared to the dispersion relation. In quadrant Q 2 , a similar kinematic reasoning leads to an increase of the angular frequency ω. This reasoning could be extended to quadrants Q 3 and Q 4 .

Sweeping by Horizontal Velocity
In this second case of uniform background flow, we consider a horizontal convecting velocity field c = (c x , 0, 0). The wave propagation is therefore convected to the negative x direction, and the initial vertical axisymmetry of the flow is broken. As a result, the dispersion relation is observed in Figure 5 to be modified proportionally to c · k = c x k x .
Unlike the case of vertical convecting velocity, the sweeping effect of horizontal convection velocity cannot be computed by an exact analytical solution in the (ω, θ, k) plane. It can only be estimated within an interval between the minimum and maximum of the analytical dispersion relation. This is because the azimuthal angle ϕ is discarded when computing the energy density in Equation (5), due to the fact that the projection of the wavevector k along the x axis is k x = k cos θ sin ϕ (see Figure 1a). At fixed radial angle θ and fixed wavenumber k, all angles ϕ = 0 · · · 2π are selected and thus k x may vary from 0 to k h . The region within which the frequency ω is found therefore lies between two boundaries: The modified Saint Andrew's cross pattern advected by a horizontal velocity c x = −10 −2 is shown in Figure 5a for the same physical and numerical parameters as in Section 2.4. In physical space, Figure 5a shows that the wave packet propagation pattern is different from that of the vertical velocity case (Figure 5b): the flow being axisymmetric along the axis of gravity, a sweeping velocity in the vertical direction preserves this symmetry, but horizontal sweeping breaks it. In frequency space, Figure 5b shows the distribution of energy density E(θ, ω): zones of maximum energy indicate that the dispersion relation is modified consistently with the above-mentioned enveloping region defined in Equation (11), with maximum modification with respect to the dispersion relation when 56 < k ≤ k max = 60. In addition, we note that the intersection between the forcing frequency horizontal line and the enveloping region, the energy density E(θ, ω, 56 < k ≤ 60) blows up, indicating maximal response of the linear system dynamics. Note that Figure 5b is in agreement with Figure 1b: in quadrant Q 1 , c is in same direction compared to the horizontal component of the phase velocity v φ for the wavevector k, and then the frequency ω increase compared to the dispersion relation by a factor that depends on the alignment of v φ with c, i.e., the azimuthal angle ϕ. The same reasoning applied to the wavevector k leads to a decrease in the frequency ω compared to the dispersion relation ω r . Therefore, in quadrant Q 1 , the frequency of the waves surrounds the dispersion relation ω r . The same kinematic reasoning can be done for quadrants Q 2 , Q 3 and Q 4 .

Doppler Shift Due to the Motion of the Wave Source
In Section 3.1, we consider how the wave frequency could be modified by the presence of a convective velocity. We now investigate a reversed configuration in which the wave production mechanism moves at a constant velocity, and in which the modification of the dispersion curve could be different. This new configuration is simulated using a forcing F b modified so that it is spatially phase-shifted, which, in practice, is done in Fourier space. In physical space, the forcing is thus advected at constant vertical velocity c z . The new forcing is therefore: with ω f = 0.3 and we choose c z = 6.28 × 10 −3 .
In Figure 6a, in the physical vertical plane, the local forcing can be seen to move in the positive z direction in physical space for the same physical and numerical parameter as in Section 2.4. The wave seems to be similar to the case of the sweeping effect with vertical velocity of Figure 4a, with an opposite direction of propagation.
In the Fourier domain, the dispersion curve is not modified, but only the forcing frequency ω f is modified by addition of the correction c · k = c z k sin θ. It is modified similarly to the sweeping effect as the new forcing frequency (dash-dotted yellow curve in Figure 6b) is equal to Therefore, unlike the sweeping effect, the Doppler effect does not modify the frequency of internal gravity waves and the original dispersion relation is preserved. This is shown in Figure 6b, where the dispersion relation Equation (2) is not modified (red dashed curve), whereas the forcing frequency is modified according to Equation (13) for a wavenumber amplitude between 56 < k < k max = 60. At the crossing between the two curves indicating the modified frequency of forcing and the dispersion relation, the concentration of energy density E(θ, ω) is maximal at (θ, ω) = (−1.5, 0.1) and (−0.9, 0.6), and their symmetries with respect to the (0, 0) point are at (0.9, −0.6) and (1.5, −0.1)

Sweeping by a Sheared Convective Mean Flow
The propagation of waves convected by a large-scale non uniform flow is considered here. In a stratified flow, a large scale shear flow can often dominate the overall structure of the flow, and evolves slowly in time. It is called the Vertically Sheared Horizontal Flow (VSHF). This shear flow is characterized by vertical variation only, without vertical velocity. This strong vertical shearing depends only on a vertical wavevector (i.e., k h = 0), so that the velocity u shear (k h = 0, k z ) is parallel to the horizontal plane and the dispersion relation gives a frequency ω = 0 (see [30]). This VSHF mode is not considered as a wave and the perturbation density distribution b in this case is uniformly zero. Therefore, these large-scale non-uniform flows are often the dominant sources of sweeping of internal gravity waves. This mode has generally been observed in numerical simulations, either in forced turbulence [30][31][32] or in decaying turbulence [33], as well as in statistical approaches such as EDQNM models [21] or the statistical state dynamics [12]. Moreover, VSHF appear in strongly stratified turbulence (low Froude number and high buoyancy Reynolds number) and not in weak stratification (high Froude number) [34]. In the context of geophysical fluid dynamics, the VSHF is an ageostrophic horizontal wind and can help to understand the emergence and maintenance of some turbulent jets, such as banded winds of Jupiter or equatorial deep jets of ocean by an interaction between wave and mean flow [35]. Nevertheless, the formation of VSHF is not really understood in many ways. Several theories try to explain the mechanisms of VSHF's formation, such as resonant interactions among gravity wave [30,32], rapid distortion theory [36] and recently by a linearized version of statistical state dynamics [12]. It appears that the VSHF is fed by an interaction with small scales of turbulence, so that a minimum rate of turbulence is required. To quantify the effect of the VSHF mode on the dispersion relation and on the internal gravity waves propagation, our method is to proceed in two stages.
In the first stage, the VSHF flow is obtained from a direct numerical simulation of the linearized Navier-Stokes-Boussinesq Equation (1). The VSHF mode is directly forced using a stochastic forcing of velocity explained in [28]. The force is only imposed on the velocity component F u (k h = 0, k z ) with k z = ±1 and an injected power P =< F u · u >= 10 −6 ; as a result, the forcing only operates on the large-scales of the flow. The obtained VSHF is a solution of the linearized Equation (1). Figure 7a shows the distribution in the (x, z) vertical plane of physical space of the horizontal component u shear,x of velocity u shear = (u shear,x , u shear,y , 0).
In the second stage, the amplitude of the VSHF flow is chosen to produce a significant sweeping effect although without totally drowning the original dispersion relation. This ad hoc adjustment permits to observe the qualitative effects of VSHF sweeping on internal waves propagation, but should eventually be related to physical observation or experimental data for better quantification. In the following, the shear velocity u shear is increased by 2.5 to obtain a maximum for the horizontal shear velocity norm of max(|u shear |) = 0.01.
We then consider the flow containing the obtained VSHF mode u shear , and its linear dynamics is numerically solved using the following linearized Navier-Stokes-Boussinesq equations: where ω shear = ∇ × u shear is the mean vorticity induced by the VSHF. Equation (14) is solved numerically with a resolution n grid = 128 grid points in each space direction, ∆t = 0.25, ν = 10 −7 , ω f = 0.5, and by using a point forcing as in Section 2.4. The simulation is started from fluid at rest, and, from this initial time on, statistics are accumulated, so that the analyzed dataset is composed of 3000 velocity-buoyancy fields equally spaced in time over 119.4T N where T N = 2π/N is the Brunt-Väisälä period.
In Figure 7b, the cross-shaped propagation pattern is distorted and does not follow a straight propagation line. The direction of propagation of the wave changes along with the advection velocity it encounters locally. This shows that the frequency of the waves is modified along its propagation resulting in a different angle of iso-velocity lines.
In Figure 8, the dispersion relation is clearly modified due to the VSHF mode both in the small wave number range 16 < k ≤ 20 ( Figure 8a) and in the large wavenumber range 56 < k ≤ 60 ( Figure 8b). On these plots, we also indicate an estimation of the sweeping effect amplitude based on the product of the amplitude of the maximum shear velocity vector with the wavevector k, in the region described by the following two boundaries: where the convective velocity is c x = rms(u shear,x ) and c y = rms(u shear,y ). According to Figure 8, a concentration of energy can be found at the intersection of the wave frequency and the modified dispersion relation, at both large-scales (16 < k < 20) and small scales (56 < k < 60). To conclude, there is an amplification of forcing within the region defined by Equation (15): the larger is the wavenumber amplitude k, the larger is the sweeping effect. Moreover, the wave frequency, estimated by the sweeping effect only, seems to be in good agreement with the wave frequency observed in the simulation, notwithstanding the neglect of the effects of the term ω × u shear + ω shear × u in the modification of the dispersion relation. This suggests that the sweeping effect is the main source of modification of wave frequency. Therefore, to estimate the wave frequency in a fully turbulent flow, one might simply estimate the value of the advective velocity responsible of the sweeping effect.

Four-Dimensional Filtering
As mentioned in the Introduction, the classical method of decomposing stratified flows between "waves" and "vortices" [17] is a spatial decomposition. This method assumes that vortices are only on horizontal plane (in the toroidal direction), whereas waves are only in the poloidal direction (see Figure A2). This is true for a linear problem, where only the dispersion relation frequency is excited in the poloidal direction. Nevertheless, for a complete non-linear problem, the poloidal component contains some eddy motion, which has a different frequency from the dispersion relation frequency. Therefore, according to this method, they are considered as "waves" for all frequencies of the poloidal component, even those outside the dispersion relation. In that case, the explicit time dependence of velocity modes must be taken into account to discriminate waves from turbulent motion. In this case, a part of the poloidal component is considered to be an eddy part, and they are now able to vertically overturn the density field as opposed to the toroidal component. For that, a first solution is proposed in Section 3. The method of detection of waves exposed and illustrated in Section 3 is based on reduced statistics to obtain the spectral distribution of energy in θ and ω. The cost of this method is limited, and it provides a first indication of the amount of wavy motion in the flow.
If one wishes to proceed to a more advanced decomposition, the solution is therefore to compute the complete space-time transforms of velocity-buoyancy fields, and to rely on this detailed information for estimating the space-time coherence of motion and thus detect the presence of wave packets. This is what we call 4D analysis, because four-dimensional transformations are required: three dimensions for space and one dimension for time. The main idea of our decomposition is to filter in Fourier space what we call "waves" by selecting the Fourier modes corresponding to the expected modified dispersion relation of internal gravity waves. The remaining part of the flow field shall be considered as "vortices", or more broadly as eddy turbulence, and from here on we call it "eddy" to distinguish with the previous classical wave-vortex decomposition by Riley et al. [17]. Note that the dispersion relation is modified by a unique sweeping velocity in the laminar case (see previous Section 3.1), whereas, in the turbulent case, there is a continuum of sweeping velocities that we assume to vary between −c and +c. It is not easy to estimate a relevant velocity c = (c x , c y , c z ) which is representative of the sweeping effect in the three spatial directions. We detail this estimation in Section 4.4. The algorithm for the new wave-eddy decomposition is detailed hereafter in five steps (see also schematic in Figure A3): 1. Direct Numerical Simulation (DNS) provides time-dependent velocity-buoyancy spectral coefficients in 3D Fourier space in terms of the wavevector k. In the following, we denote by f either the velocity or buoyancy coefficient, as f (k, t). 2. We compute the 1D time-Fourier transform as FFT( f (k, t)) = f (k, ω). 3. For each point (k x , k y , k z , ω) in the 4D space, the "waves" contribution to the flow field is denoted f w (k, ω), and is obtained by filtering the Fourier modes matching the modified dispersion relation, i.e., lying in the region defined by: where the range of advective velocities is characterized by a single c = (c x , c y , c z ) is obtained from DNS. The "eddy" part is denoted f e (k, ω) and obtained as the difference between f and f w . For the horizontal velocity field, the VSHF mode (ω 0 and k h = 0) is removed from the wave component. It corresponds to the point (θ = ±π/2, ω 0) in the (θ, ω) plane. 4. We compute the 1D inverse Fourier transforms FFT −1 ( f w (k, ω)) = f w (k, t) and FFT −1 ( f e (k, ω)) = f e (k, t) to return to physical time for the wave and eddy coefficients. 5. Compute the 3D inverse Fourier transform in physical space, which permits recovering the wave and eddy fields in physical space and to visualize them.
Note that, in Step 3 of the filtering, the VSHF part of the horizontal velocity of the wave component with an angular frequency ω 0 and with a horizontal wavenumber k h = 0 is removed. It does not appear either in the eddy component but it appears in the full velocity field. This is done because the VSHF would dominate the horizontal velocity field of the wave. Of course, this algorithm is computationally demanding in terms of both CPU time and memory storage.

Application to the Sweeping by Horizontal Velocity of Saint-Andrews cross Wave Propagation
In this section, the previous technique presented to separate the wave and the eddy components of the flow in Section 4.1 is applied on the same numerical simulation as in Section 3.1.2. The numerical parameters are also the same as in Section 3.1.2. Therefore, the chosen convective velocity is c = (10 −2 , 0, 0). As the convective velocity is in only one direction, one could decide to refine the filtering in Step 3 of Section 4.1 to take only one direction into account. However, for simplicity, the filtering technique is kept as if the convective flow were in both positive and negative x direction. According to Step 3 of the method, we computed the vertical velocity u z (k, ω) in four-dimensional space (k, ω). Since the 4D filtering is defined for every (k, ω), the velocity is either a wave part u z = u w z or an eddy part u z = u e z , depending on the value of (k, ω). Then, the three vertical densities of the energy are defined (denoted with to distinguish them with the previous energy given in Equation (5) and associated to u z , u w z and u e z : E ,e,w z (θ, ω) = ∑ θ(k)∈I θ ,0≤|k|≤k max 1 2 |u e,w z (k, ω)| 2 (18) with this definition, E z (θ, ω) = E ,e z (θ, ω) + E ,w z (θ, ω). In Figure 9, one can see E z for the full vertical velocity field (Figure 9a), the wave part E ,w z (Figure 9b) and the eddy part E ,e z (Figure 9c). In Figure 9c, we see that there is no energy along the dispersion relation (indicated with a red dashed line). This is the imprint of the discretization error of the filtering technique. Furthermore, the density of vertical energy is quite low in Figure 9c when, at the same location in (θ, ω) plane in Figure 9b, the density of vertical energy is high. Indeed, most of the energy in this area comes from the waves that have been modified in angular frequency by the sweeping effect (for large wavenumber k x , i.e., small scale). It is the manifestation of the density of energy sharing the same coordinate in the (θ,ω) plane, but for a different wavevector k and consequently a different sweeping effect.
Compared to the total energy in Figure 9b using the first separation method, our 4D analysis appears more precise-i.e., without vertical line-whereas no treatment such as Hann windowing is applied to the data. The main reason for this is due to the difference in the two algorithms for obtaining the energy density: The first algorithm gathers the flow component in a θ ∈ I θ domain before doing the time Fourier transform. It is only after applying the Fourier transform that the energy density is computed. To understand the link between the two methods, by analogy with the total energy defined previously in Equation (5), we define the vertical energy by: The second technique, which is more computationally expensive, starts with the time Fourier transform of the flow component, then summing them in the θ ∈ I θ domain and finally computing the energy density. As shown in [23], the difference between the two formulations of vertical density of energy is: with k a different wavevector from k sharing the same domain I θ . The second term on the right-hand side corresponds to the error done by the first technique against the second one. The error corresponds to the crossing energy sharing the same domain I θ . This difference can be summed up in one sentence: the energy of the sum of the flow component in a domain I θ is different from the sum of the energy of the flow component in the same domain I θ . The final result ends up being slightly different for laminar cases. For a turbulent case, we assume that the correlation term < u(k, t)u(k , t) > is almost zero on average. After analysis in the domain (k, ω), we come back to physical space (x, t), where the vertical velocity u z (x, t) is decomposed into wave and eddy parts u z (x, t) = u e z (x, t) + u w z (x, t). Figure 10 shows u e z (x, t) and u w z (x, t). In the eddy component (Figure 10b), only the Dirac in space is visible, whereas for the wave component (Figure 10a) only the wave propagation is visible. By comparing the sum of Figure 10a,b, the same kind of flow is obtained as in Figure 5a. Therefore, we can safely assume that, for the case of simple, laminar, convective flow, this filtering technique works rather well. In the following, this analysis is applied to a DNS at moderately low Froude number and at low buoyancy Reynolds number.

DNS and Numerical Parameters
In contrast with the previous characterization of Section 3, in which we analyze the effect of large-scale mean field on the propagation of internal gravity waves produced by harmonic monochromatic forcing, and its trace on the statistics in spectral space, we consider here stratified turbulence obtained by Direct Numerical Simulations of the full nonlinear Navier-Stokes-Boussinesq equations. The corresponding field is then subjected to analysis by four-dimensional Fourier transform (FFT4D). As proposed by Maffioli [28], the flow is produced by forcing the momentum conservation equation, until statistical stationarity is reached. As presented in Section 2.2, the force F u is applied at very large-scale in a cylindrical shell of radius k h = 4 and vertical extent 1 ≤ k v ≤ 3, thus the forcing wavenumber is approximately k f orc 5 and the angles of forced wavevectors is in the range 0.24 rad < |θ| < 0.65 rad. This means that the VSHF is not directly fed by the forcing, so that its growth is limited, even if it eventually becomes energetically dominant in the flow. The injected power is P =< F u · u >= 0.01. The resolution is n grid = 256 points in each direction, viscosity is ν = 1/2500, the Brunt-Väisälä frequency is N = 5 and the time step is ∆t = 0.005.
The simulation is run until a statistically stationary state for total kinetic energy is reached. In this state, the Froude number is Fr = ε/Nu 2 h = 0.07, where ε is the dissipation of kinetic energy and u h = (K h − K shear ) 1/2 is the rms of horizontal velocity component with K h the horizontal kinetic energy and K shear the shear mode energy (i.e., energy of all the k h = 0 Fourier modes). The Kolmogorov length scale is η = (ν 3 /ε) 0.25 = 9 × 10 −3 and the Kolmogorov time scale is τ η = (ν/ε) 0.5 = 0.2. The buoyancy Reynolds number is Re b = ε/νN 2 = 1. The integral length scale L = 0.4 is rather large, resulting in a flow dominated by waves and explaining why the Taylor-length-scale-based Reynolds number is high as well Re λ = 95. The Ozmidov scale is k 0 = (N 3 /ε) 1/2 ≈ 100 ≥ k max , so that wave motion is expected to be mixed with turbulent motion at all flow scales.
The five-step algorithm for spatiotemporal decomposition is applied on velocity and buoyancy using 1000 snapshots of these fields separated by a time step of ∆t = 10∆t = 0.05. The total time span is therefore T = 39.8T N where T N = 2π/N is the Brunt-Väisälä period. Due to this discretization of frequency with a step ∆ω = 2π/T, the angular frequency ω is determined to within ∆ω accuracy, or equivalently the angular frequency cut-off is ∆ω. The filtered range described in Equation (16) is therefore also defined within a 2∆ω accuracy.

Estimation of the Sweeping Velocity
As explained in Section 3.3, the VSHF can sometimes dominate the overall structure of the flow. To check this in our numerical simulation, the kinetic energy of the VSHF flow K shear based on horizontal velocity u shear = u(k h = 0, k z ), and the total kinetic energy K tot are obtained as: Figure 11a shows the time evolution of K shear and K tot from the beginning of the simulation, where initial conditions with zero velocity and density fluctuations were used. The statistically steady state is reached at t 30, and there is a slow growth of energy at the end of the simulation, due to the growth of VSHF. The kinetic energy of the VSHF amounts to about a third of the total kinetic energy, so that the VSHF can be considered as the main source of modification of wave frequency by sweeping. To obtain the sweeping velocity related to the shear, one computes a rms shear velocity in the x and y directions (u shear,x , u shear,y ) from the shear kinetic energy K shear : u shear,x = u shear,y = K shear ≈ 0.14 .
Therefore, the sweeping velocity is estimated as c = (u shear,x , u shear,y , 0) and in the numerical simulation, we use c = (0.14, 0.14, 0) from Equation (23). Choosing a different sweeping velocity would inevitably change the result of the filtering, so that the choice of the sweeping velocity is probably the main source of error of this technique. Here, we think that most of the sweeping effect comes from the VSHF mode, but this is suggestive.

Time-Windowing for the Time Fourier Transform
The choice of parameters for the time span over which Fourier transforms are done can be of importance, since a too narrow time-window or too sparse time points can lead to significant error in angular frequency estimates. Therefore, after transforming the velocity and buoyancy fields back to physical time t domain, we compute the time-evolution of the wave and the eddy contributions to vertical kinetic energy: These energies are plotted in Figure 11b. We first observe that the energy in the wave part is about six times larger than the energy in the eddy part. Second, we see that for most of the time domain, away from the boundaries t = 50 or t = 100, the energies evolution only consists of fluctuations around a mean value. The appearance of edge zones in the evolution is due to the filtering technique which uses frequency-separated components that break up the periodic structure of the velocity field and result in a transition period at the edge of the time domain. However, these edge zones are limited in extent, so that we use no windowing technique such as Hann windowing, as often done in, e.g., experimental data processing. Furthermore, after a windowing technique, to recover the amplitude or the density of energy correctly, the signal needs to be modified by a different factor which might lead to inaccuracies later on in the energy spectrum and the visualization of the flow components. In the following, the presented statistics are obtained using only the part of the time window away from the edge zones.

Energy Distribution in the (θ, ω) Plane
The filtering technique is applied on the four variables f = u x , u y , u z , and b obtained from the DNS, by the five-step algorithm method explained in Section 4.1. The representation of the results requires a reduction of the amount of information contained in the 4D-fields depending on (k, ω). Therefore, the fields are averaged over wavenumber k and azimuthal angle ϕ to retain only the two dimensional dependence on (θ, ω). Thus, by analogy with the vertical density of energy defined in Equation (17), we define the averaged horizontal and potential energies by: These quantities are not the same as the energy density defined in Section 2.3. The 4D-transform technique is more expensive, but more accurate in terms of determination of the wave contents of the flow than the method of Section 2.3. Figure 12 shows the results obtained when applying the 4D-decomposition technique to the results of our simulation. Figures 12 first column show the distributions of horizontal, vertical and potential energies in the (θ, ω) plane. Figure 12 (Column 2) shows the distributions of their wave contents as obtained after the separating procedure and Figure 12 (Column 3) shows the distribution of the eddy contents of the same energies. Figure 12 (Column 1)-for the complete horizontal, vertical and buoyancy energies E h , E z and E b -shows that, after the spatiotemporal decomposition, the energy is concentrated on the dispersion relation curve of the waves, as well as away from the dispersion relation, for high angle θ and for a large range of angular frequency ω. We observe a concentration of horizontal and buoyancy energies along wavevectors ks almost vertical (θ close to ±π/2), which can be attributed to horizontal structures. In addition, Figure 12a exhibits a concentration of horizontal energy density along the line ω 0, which indicates motion unrelated to the waves, thus could be attributed to the eddy contents of the flow. The structure of horizontally-stretched pancake-like eddies comes to mind, as observed in previous studies (e.g., the experimental study by Fincham et al. [37]). (17), (25) and (26) of horizontal and vertical kinetic energies and of potential energy, we define the same energies for the wave and eddy contents of the flow by using the analogy with vertical density of energy Equation (18):

Similar to the definitions in Equations
• For horizontal and potential energy density, we define the wave part of the energies, respectively noted E ,w h (θ, ω), E ,w b (θ, ω), based on u w x ,u w y , and b w .

•
Similarly, we define the eddy part of energy, with energies noted E ,e h (θ, ω), E ,e b (θ, ω), based on u e x , u e y , b e .
These quantities are plotted in Figure 12 (Columns 2 and 3). By comparing the wave part ( * w ) and eddy part ( * e ), our decomposition separates clearly the wave energy, corresponding to the condensation of the energy around the dispersion relation curve (Figure 12, Column 2), from the energy associated to more general turbulent motion. This other motion is either the eddy part of the flow, with a concentration of energy around ω 0 in Figure 12c, or horizontal structure which concentrates energy along vertical wave number regions (Figure 12c,i). The energy of the waves appears to be concentrated around the forcing angle 0.24rad < |θ| < 0.65rad, and decreases only for E ,w z (θ, ω) where θ ±π/2. The overall energy density level in Figure 12f is clearly lower than that of the other figures, since it corresponds to the vertical kinetic energy of the non-wavy and non-VSHF part of the flow.
Moreover, some energy appears in the neighborhood of the relation dispersion curve (light blue part in Figure 12, Column 2) probably due to wave numbers that undergo a strong sweeping effect due to their small scale. For example, the relative sweeping amplitude is larger for a close-to-horizontal wavevector with large wavenumber-such that k x k z 1 so that θ 0-than for a wavevector at smaller wavenumber such that k x ∼ 1 k z ∼ 0 for which θ 0. The strong sweeping effect can also be seen in Figure 12 (Column 3) where the edge of the density of energy along the dispersion relation is irregular. Indeed, at some angle θ, a small value of k x ∼ 0 and k y ∼ 0 might not exist in the spectral discretization, therefore resulting in a minimum sweeping effect already large as it is multiplied by a large value (c x = c y = 0.14). For another angle θ, there might exist smaller values of k x and k y , which result in a minimum sweeping effect smaller as well. This effect does not exist in Figure 9b because the sweeping velocity is small c = (10 −2 , 0, 0).

Energy Spectra
The wave-eddy decomposition proposed above can be used to obtain several statistics of different nature for the separate wave and eddy contents.
We focus here on the energy spectrum at scale k for the wave and eddy parts, which can be computed for horizontal, vertical and potential energies, respectively, as The energy spectra are plotted in Figure 13. At large-scale, the eddy part is dominant for horizontal energy (Figure 13a), whereas the wave part is dominant for both vertical (Figure 13b) and potential energies (Figure 13c). This is consistent with results in the previous section, wherein the accumulation of energy density along the dispersion relation curve is more important in vertical and potential energy, whereas the neighborhood of frequency ω ∼ 0 accumulates more eddy energy. However, at smaller scales, there seems to be an equipartition of eddy and wave energy in the horizontal energy spectrum and the buoyancy energy spectrum. Moreover, Figure 13b for the vertical kinetic energy spectrum shows that wave energy is dominant even at small scales.
Since the Reynolds number of our simulation is relatively low, one cannot expect to observe clear inertial scalings. We note however that the wave and eddy energy spectra seem to exhibit scalings respectively close to k −3 and k −5/3 for both the horizontal and the vertical kinetic energy spectra. Nonetheless, this is still suggestive and more remains to be done to understand what the correct scaling should be.
Finally, Figure 13 shows several spectra in the time interval t ∈ [50; 100], where the flow is statistically steady in terms of total integrated kinetic energy. It is not possible to distinguish between the different curves because they are on top of each other, but by looking at the width of each line, one can understand how much the energy spectra change through time. The different curves for each kind of spectrum show that there is some variability in the spectra with time, especially in the largest scale in the infrared spectrum range. However, the spectra are still remarkably stable on average over the chosen period.

Physical Space Distribution of Wave and Eddy Fields
The final use of our decomposition for the current DNS results is the physical space distribution of the field. Upon computing the 3D spatial inverse Fourier transform of the separated wave and eddy spectral coefficients, one obtains the velocity and buoyancy fields u DNS x , u DNS y , u DNS z , and b DNS decomposed in terms of their wave and eddy parts u e,w x , u e,w y , u e,w z , and b e,w in physical space. Figure 14, shows their distribution in the vertical (x,z) plane at mid y-span. Therefore, the wave and eddy parts of the velocity field and of the buoyancy field are recovered. In the first column of Figure 14, the full field is plotted as for classical 3D DNS results. In the second column of Figure 14, the distribution of the wave contents of the velocity field is represented, and, in the third column of Figure 14, that of the eddy part is shown. The time evolution of these fields is informative of their proper kinematics. Therefore, the animation of these three decompositions for the three variables in the time interval t ∈ [60; 90] is available in the Supplementary Materials (see Videoa S1-S3). We recall that the shear mode is not taken into account either in the wave part or in the vortex part.
In Figure 14 (Rows 2 and 3), we represent the vertical velocity field u z (x, z) and buoyancy field b in (x, z) plane and decomposed into a wave and an eddy parts. It appears that the wave part is at a bigger scale than the eddy part of these two fields. This underlines that, in the Energy spectrum of Figure 13b,c, the wave energy dominates the eddy energy at larger scale. In addition, their time evolution shows that the wave part of u z (x, z) and b oscillates like a wave while it is not the case in the eddy part (Videos S2 and S3). On the other hand, the decomposition of the horizontal velocity u x (x, z), which is illustrated in Figure 14 (Row 1), shows clearly the overall horizontal layering either for wave and eddy parts. These layers seem to be at a smaller vertical scale for the wave part than for the eddy part of the u x field.
These observations match with the Energy spectrum of Figure 13a, wherein the eddy energy dominates the wave energy at large scale. Again, the time evolution shows that the wave part of u x (x, z) oscillates as a wave while it is not the case in the eddy part. The layers in the wave part u w x appear surprising, even though, intuitively, the layering is associated with the eddy part u e x . Indeed, u w z contains only the vertical component of the poloidal part of velocity (see Appendix B and Figure A2) which is associated to linear wave, whereas u w x contains a mix of horizontal components of the poloidal and toroidal part of the velocity, which are associated with the linear wave and the horizontal vortex respectively (see Appendix B and Figure A2). The direct analysis on the horizontal velocity appears as a limitation of our method because we call wave part u w 1 any mode of u 1 which has a frequency and a wavevector included in the dispersion relation without distinguishing the poloidal and toroidal parts in u x . It would be preferable to use our method directly on the poloidal and toroidal components separately. This limitation occurs in experimental work where the evolution of a three-dimensional field is generally not possible.
From a qualitative point of view, this new decomposition should work when the domain of the waves and the domain of the eddies are disjoint in (k, ω) space. One expects the angular frequency of the eddies to be small (i.e., ω 0, see Figure 13c) whereas the angular frequency of the waves are engulfed between the positive and negative sweeping effect (ω r − |c · k| < ω < ω r + |c · k|). Therefore, if the sweeping effect due to the advective velocity c is strong enough, then some waves with high amplitude of wavevector have a strong deviation of frequency by the sweeping effect reaching a zero frequency (ω = N cos θ + c · k 0) and therefore sharing the same frequency as the majority of the eddies. This happens for the higher wave number k η = 1/η based on the Kolmogorov scale η. From this observation, one could compute a dimensionless number: the ratio of the Brunt-Väisälä frequency N over the sweeping effect k η c: N/(k η c). If this ratio is higher than one, the domains of the waves and eddies are well separated in (k, ω) space: ω 0 for the eddies and ω > 0 for the waves. If the ratio is lower than one, the domains of the waves and eddies overlap in (k, ω) space. In this study, this dimensionless number would be close to 0.2 with c 0.2, N = 5 and k η 110. This might explain why the approximation of the horizontal wave component is only approximative, as most of the horizontal energy comes from the motion of the eddies. The ratio N/(k η c) ≥ 1 could be obtained by changing N or decreasing c for the same resolution. Nevertheless, apart from these limitations on the horizontal velocity field, the analysis of vertical velocity field u z (x, z) and buoyancy field b provide valuable information about the overturn of density.
We recall that in the standard wave-vortex decomposition [17] all the buoyancy field b is attributed to wave motion. This is not the case with the method we propose, and we can classify the density turnover in the eddy part. We illustrate this distinction in Figure 15 with a zoom in of the temporal evolution of the buoyancy field in the vertical (x, z) plane, for one overturning event. We observe that the eddy part of b (Figure 15 first line) shows an overturn of density [38] by looking at the black iso-density lines: when the iso-density lines are on top of one another (as in Figure 15b), heavy fluid is on top of lighter fluid. Consequently, mixing is occurring as heavier fluid tends to flow below lighter fluid. On the contrary, in Figure 15 (Row 2), which corresponds to the wave part of the buoyancy field, no similar topological arrangement is seen. However, the amplitude of iso-density lines in the middle of Figure 15d is quite high, whereas, in Figure 15e,f, it is lower. This could suggest that the waves are interacting and are breaking into eddies, but this is still speculative. Other events of this type are visible in the animation and no density overturn is visible on the wave part of b. Nonetheless, this shows that the present wave-eddy decomposition successfully separates a turbulent event, the overturning of density, from the wave component.
By performing this decomposition in physical space, it is then possible to extend the analysis to other statistics in physical space (structure function, two-point correlation, etc.) on each part separately in order to obtain information on the role of each part and their interaction.

Conclusions and Perspectives
In this work, we propose numerical simulations of Navier-Stokes-Boussinesq for stably stratified flows generated by specific forcings, in which we use two spatiotemporal techniques for separating the wave motion from the eddy motion. The statistical characterization is mostly based on the computation of the energy density in the (θ, ω) plane.
The first part presents a spatiotemporal analysis used in simple laminar cases (homogeneous advective flow and shear flow) in order to quantify the modification by a large-scale mean flow of the propagation of the waves, in terms of wave frequency and the observed dispersion relation in the simulations. It appears that the main physical reason for the modification of the wave frequency is the sweeping effect. This effect corresponds to the advection of the wave packets by a large-scale flow, in particular by vertically sheared horizontal flow (VSHF).
The second part presents a new spatiotemporal decomposition, more precise but more numerically expensive than the first method: the partitioning of waves and eddies is done by separating the angular frequency of the waves from the rest of the turbulence. This separation takes into account the sweeping effect due to the shear part of the flow (VSHF), which is supposed to consist of a large-scale advecting velocity. From this separation, energy statistics and visualization of the velocity and buoyancy field can be done for the wave and the eddy parts of the flow.
Such a decomposition is more sophisticated than the decomposition proposed by Riley [17] in the physical space only. The latter decomposition is solely based on the physical domain and lacks the temporal dependence of the wave component to be fully effective. Therefore, the spatiotemporal technique proposed here could be more effective, especially if it is applied in complement to the spatial decomposition proposed by Riley. From there, it is possible to use advanced statistical characterization in physical space (structure function, two-point correlation, etc.) to obtain further information on the role of the different parts and their interaction. For example, our decomposition is able to distinguish a turnover of the density field associated with the eddy part from an oscillation of the density field associated with the wave part. However, it still needs to be further studied. In our present work, observations of the sweeping effect done on laminar cases are extended to the turbulent case. However, the analogy between the turbulent and laminar cases may be limited. A better understanding of the sweeping effect would therefore be needed for turbulent cases. The other main limitation of this decomposition is its dependency on the choice of the advective velocity responsible for the sweeping effect. We assume that the VSHF is the main source of sweeping effect, which seems to be a good approximation but not a completely universal one. As the result of the decomposition depends significantly on the chosen advective velocity, a better evaluation of this velocity is of importance for the accuracy of our decomposition. This technique is valid when the domain of the waves and the domain of the eddies are well separated in (k, ω) space. To quantify this separation, a dimensionless number has been introduced in Section 4.8, the ratio of the Brunt-Väisälä frequency N over the sweeping effect k η c: N/(k η c). Even though N/(k η c) 0.2 in this study, which means that the sweeping effect is strong against the stratification, waves and eddies have been well separated for the buoyancy and vertical velocity field. Therefore, further work remains to be done to understand what the limits of this new algorithm are against the dimensionless number N/(k η c).
We have shown the applicability of our method on an instance of turbulence where the wave regime dominates the eddy regime, at large Ozmidov scale and Re b 1. It will be necessary to extend our decomposition to a wide range of parameters in terms of Froude number (lower Fr) and buoyancy Reynolds number (higher Re b ). It will also be necessary to increase the Ozmidov scale and retaining enough wave contents. For this, numerical simulations are required where both waves and eddies dominate a different domain of the energy spectrum, but they would also be more computationally expensive.
If this new decomposition was successfully further developed, many new physical applications could be proposed: computing the ratio of energy between waves and eddies, or the different energy spectra against horizontal and vertical wavenumber. Another possibility would be to compute energy transfers between waves, eddies and themselves using similar techniques done in the case of magnetohydrodynamics [39]. Moreover, this spatiotemporal decomposition method could be applied to simulations of rotating turbulence for which experimental studies have already used a similar 4D Fourier technique [8], but more interestingly to mixed inertio-gravity waves in rotating and stably stratified flows.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4433/11/4/420/ s1, Video S1: Full, wave and eddy horizontal velocity component u x in the (x, z) plane in the middle of the y interval. Video S2: Full, wave and eddy vertical velocity component in the (x, z) plane in the middle of the y interval. Video S3: Full, wave and eddy buoyancy field with superimposed iso-density lines in black in the (x, z) plane in the middle of the y interval. Funding: This research was funded by ANR DisET grant number ANR-17-CE30-0003. This work was granted access to the HPC resources of IDRIS under the allocation A0062A02206 made by GENCI and HPC resources of the FLMSN, partner of EQUIPEX EQUIP@MESO. convergence of the velocity fields. On the contrary, after Section 4.3, the convergence of the kinetic energy is reached, which reduces the advantage of a windowing technique. Therefore, to stay coherent inside the entire article, we do not use any windowing technique.
In Figure A1, there are still some vertical lines that can be seen. This is due to the viscosity and to the discretization error of the forcing (the forcing is not exactly a Dirac).

Appendix B. Analytical Solution for a Local Oscillating Buoyancy Forcing
In Section 3, we observe numerically the flow obtained by the harmonic forcing (Equation (6)) of buoyancy: in physical (x, z) plane, the flow exhibits a Saint Andrew's cross pattern, and the energy concentration map in (θ, ω) plane exhibits the complete dispersion relation curve, as if waves at all frequencies were present. We explain here theoretically how it is possible to obtain the full dispersion relation from a local forcing at a constant frequency from the linearized Navier-Stokes-Boussinesq equations. We start with the linearized inviscid Navier-Stokes equation in the Boussinesq approximation: Taking the divergence of Equation (A1), one finds the Poisson equation The spatial Fourier transform of Equations (A1) and (A4) yields Combining Equations (A5) and (A6), one obtains Taking the spatial Fourier transform of Equation (A3) and multiplying by the unit vector k/k, we get: Figure A2 illustrates the polar-spherical reference frame which, in the context of turbulent flow decomposition, especially for stratified flows, was named Craya-Herring frame [20,40]. The vector e t corresponds to the toroidal component, e p corresponds to the poloidal component and e r corresponds to the radial component. In this frame of reference, the flow is incompressible (∇ · u = 0 or k ·ū = 0), which means the velocity vectorū in Fourier space is perpendicular to the wavevector k. As a result, the new velocity field in Fourier space has only two components:ū =v t e t +v p e p . Figure A2. The polar-spherical frame of reference linked to k in spectral space, also named Craya-Herring frame. Unit vectors are toroidal e t , poloidal e p , and radial e r .
Projecting Equations (A7) and (A8) on to the Craya-Herring frame, one therefore gets a new, simpler set of kinematic equations The toroidal component of Equation (A9) is no longer considered because it does not interact with the component of the other directions. The viscosity terms are also removed for the sake of simplicity. By adding a forcing on the buoyancy term, one gets the system of equations whose analytical solution we seek: Denoting ω r = N cos θ and deriving the second line using the first line, one gets ∂ 2 tb (t) + ω 2 rb (t) = ω f cos ω f t (A11)