Louis-Alexandre Couston

Daniel Lecoanet

Benjamin Favier

Michael Le Bars

Order out of chaos: slowly-reversing mean flows emerge from turbulently-generated internal waves

We demonstrate via direct numerical simulations that a periodic, oscillating mean flow spontaneously develops from turbulently-generated internal waves. We consider a minimal physical model where the fluid self-organizes in a convective layer adjacent to a stably-stratified one. Internal waves are excited by turbulent convective motions, then non-linearly interact to produce a mean flow reversing on time scales much longer than the waves' period. Our results demonstrate for the first time that the three-scale dynamics due to convection, waves, and mean flow, is generic and hence can occur in many astro/geophysical fluids. We discuss efforts to reproduce the mean flow in reduced models, where the turbulence is bypassed. We demonstrate that wave intermittency, resulting from the chaotic nature of convection, plays a key role in the mean-flow dynamics, which thus cannot be captured using only second-order statistics of the turbulent motions.

An outstanding question in fluid dynamics is whether large-scale flows can be accurately captured in reduced models that do not resolve fluid motions on small spatiotemporal scales. Reduced models are necessary in many fields of fluid mechanics, since fluid phenomena often occur on a wide range of spatial and temporal scales, preventing exploration via direct numerical simulations (DNS) of the Navier-Stokes equations. This question is of interest to, for instance, the turbulence community, which has developed closure models in Large-Eddy Simulations and Reynolds-Averaged Navier-Stokes simulations [1]; the statistical physics and geophysics communities, who aim to describe the self-organization and large-scale behavior of turbulent flows [2][3][4][5]; atmospheric and oceanographic scientists, whose goals are to provide long-time predictions of the evolution of our climate using weather-ocean models with coarse resolution [6,7].

A drastic approximation would be to assume that large-scale flows and small-scale motions are dynamically decoupled, but this is rarely the case. A number of important slow large-scale flows are controlled by rapid processes at the small scales. For instance, the 22-year cycle of solar magnetism is driven by the Sun's convective interior, which evolves on month-long or shorter timescales [8,9]; upwelling of the planetary-scale thermohaline circulation of Earth's oceans hinges on enhanced mixing events that critically depend on small-scale (∼ 100 metres) internal waves [10,11]; Jupiter's zonal jets develop from small-scale turbulence patterns due to convective heat transfers in the weather layer and deep interior [12].

The generation of a large-scale flow by turbulent fluctuations can be studied by spatial-averaging the Navier-Stokes equations. Let us consider the case of a large-scale mean flow ū in the horizontal x direction perpendicular to downward gravity. We write (u , w ) the velocity fluctuations in (x, z) directions with ẑ the upward vertical axis. In these two dimensions, the horizontal-mean of the Navier-Stokes equation in the x direction reads

∂ t ū -ν∂ zz ū = -∂ z (w u ), (1)
with ν the kinematic fluid viscosity. The right-hand side of ( 1) is minus the divergence of the Reynolds stress and is the momentum source or sink for the mean flow. In isotropic homogeneous turbulence, we do not expect the generation of a mean flow due to the lack of symmetry breaking. However, any inhomogeneity or anisotropy of the fluctuations can initiate a slowly-varying mean flow, whose fate depends on its interaction with the fluctuations [4]. The parameterization of the Reynolds stress (w u ) for unresolved scales is the key ingredient in all reduced models. Generally, a closure model expresses the Reynolds stresses in terms of the resolved variables [1].

In our case of interest, the small-scale fluctuations are oscillating disturbances of the density field called internal waves. Internal waves are ubiquituous in oceans [13], planetary atmospheres [14][15][16][17], stars [18,19], brown dwarves [20] and planetary cores [21]. In the atmosphere, internal waves actively contribute to the generation of mean equatorial winds in Earth's stratosphere, which change direction roughly every 14 months, coined the Quasi-Biennial Oscillation (QBO) [22]. Internal waves may also be involved in the generation of reversing zonal flows on Saturn [23] and Jupiter [24], they are of interest for extrasolar planetary atmospheres [25], and may influence the differential rotation of stars [26] and slow large-scale motions of Earth's magnetic field [27].

Here, we report results of the first DNS of a realistic slowly-reversing mean flow in two dimensions, and we unravel the key physics of the generation mechanism using a hierarchy of low-order models in which the Reynolds stresses are approximated. We use the horizontally-periodic self-consistent model of convectivestably-stratified dynamics of [START_REF] Couston | [END_REF]. The velocity u ¯= (u, w), temperature T , and density anomaly ρ = -αT ). Note that energy propagates upward along wave crests, so crests toward the upper left (right) correspond to retrograde (prograde) waves. (C) shows the mean flow ū(t, z). The mean flow in the convective zone corresponds to the average of stochastic plumes emitted from the bottom boundary, hence reverses on a relatively rapid, convective time scale. In the stably-stratified layer, ū results from the nonlinear interaction of internal waves and oscillates on time scales ∼ 0.1, much longer than the buoyancy period ∼ π10 -4 . Simulation details and movies are available in Supplemental Material [START_REF]See Supplemental Material for details on the models, the effect of the Prandtl number, statistics and nonlinearities[END_REF].

satisfy the Boussinesq equations

∂ t u ¯+ (u • ∇)u ¯+ ∇p = P r∇ 2 u ¯-P rRaρẑ -u ¯τ, ( 2a
)
∂ t T + (u • ∇)T = ∇ 2 T, (2b)
∇ • u = 0, (2c)
non-dimensionalized with κ (thermal diffusivity) and H (characteristic height). The fluid is thermally stratified (T t and T b imposed on the top and bottom no-slip boundaries) and exhibits a buoyancy reversal at the inversion temperature T i with T b > T i > T t (similar to water whose density maximum is at 4 • C [30]). Thus the fluid spontaneously organizes into a lower, nearly isothermal convective region, and an upper stably stratified region. P r = ν/κ and Ra = α s g∆T H 3 /(κν) are the Prandtl and global Rayleigh numbers; α s is the expansion coefficient for T > T i ; and ∆T > 0 is the difference between the dimensional bottom and inversion temperatures, such that using T i = 0 as the dimensionless reference temperature, we have T b = 1. The buoyancy reversal is obtained using the nonlinear equation of state for ρ:

ρ(T ) = -α(T )T = -T, T ≥ T i = 0, ST, T ≤ T i = 0, (3)
with S the stiffness parameter [START_REF] Couston | [END_REF]. We define the neutral buoyancy level z N B to be the height where adiabatic plumes emitted from the bottom boundary become neutrally buoyant. This corresponds to the height of the convection zone [START_REF] Couston | [END_REF]31], or equivalently, the base of the stable layer (dashed lines in figures 1A-B). The normalized domain lengths are L x = 2, L z = 1.5 in the x, z directions, which leads to an aspect ratio of the convection at statistical steady-state close to 3 for all simulations; τ = 10 2 √ 2{tanh[(z -L z + 0.15)/0.05] + 1}/2 is a z-dependent linear damping used to prevent wave reflections from the top boundary. We solve equations (2) via DNS using Dedalus [START_REF] Burns | Dedalus: A Flexible Pseudo-Spectral Framework for Solving Partial Differential Equations[END_REF] with Chebyshev polynomials (Fourier modes) in z (x) direction. DNS are run over several thermal diffusion times in order to allow the system to reach a statistical equilibrium self consistently, and obtain several reversals of the mean flow.

Figure 1 shows the main DNS results of the paper, obtained for T t = -43, T b = 1, P r = 0.2, Ra = 1.2 × 10 8 and S = 1/3, such that the convection-wave coupling is relatively strong and the interface is flexible [START_REF] Couston | [END_REF]. With z N B = 0.68, the effective Rayleigh number is Ra eff = z 3 N B Ra ≈ 4 × 10 7 . Snapshots of vertical velocity (figures 1A,B) reveal large convective updrafts and downflows below z N B , and internal waves above. If there was no mean flow in the stably stratified layer, convection would generate prograde and retrograde waves with similar amplitude. However, in figures 1A,B, the internal waves are mostly propagating in a single direction, an indication that the mean-flow (shown by the solid line) is filtering waves going in the opposite direction. The evolution of the mean flow over one thermal time scale is shown in figure 1C. The stable layer has a strong meanflow which reverses every ∼ 0.05 thermal time. Each new mean-flow phase starts near the top of the domain and descends toward the convective layer. The mean flow is driven by wave damping at critical layers and by viscous and thermal dissipation. Critical layers are ubiquitous in our DNS because convection generates a broad spectrum of waves, some of which have low phase velocities. Viscous and thermal dissipation effects are relatively strong in our DNS, so the basic mean-flow mechanism is essentially due to wave dissipation.

Previous studies of wave-mean-flow interactions have focused on momentum-deposition by internal waves of a single frequency and wavenumber [START_REF] Plumb | [END_REF]34]. In such cases, it can be shown analytically that a slowly-reversing mean flow emerges provided that there are both prograde and retrograde waves, as well as an initial disturbance. The prograde (resp. retrograde) wave provides a +x positive (resp. negative) acceleration for the meanflow through damping. Then, the competition of the two forces (whose intensity depends on the direction of the mean flow through the Doppler shift) leads to the observed long-time oscillation of ū [22].

Our results demonstrate for the first time that an oscillating mean flow can emerge from internal waves generated by turbulent motions with no control over the waves (i.e. no parameterization). Importantly, the fundamental mechanism that applies for monochromatic waves also applies for a broadband spectrum of internal waves: damping and momentum deposition is stronger for waves going in the same direction as the mean flow. This is shown in figure 1A where a strong mean flow in the positive direction strongly dissipates prograde waves, such that only retrograde waves can be visible above. The same is true in figure 1B but for the case of a negative mean flow. With a broadband spectrum of waves, whose amplitudes can vary over time due to the chaotic dynamics of convection, momentum deposition cannot be simply traced back to a handful of self-interaction terms in the Reynolds stress that would be coherent over long times. Driving of a mean flow in this context may be unexpected, but is in fact generic at sufficiently low P r: as figure 2 shows, the mean flow becomes stronger and more regular as P r decreases. This can be understood from the fact that while the forcing through wave damping is only weakly affected by decreasing P r (because waves are damped through both viscous and thermal dissipation effects), the mean flow experiences much less dissipation (it is only damped through viscosity effects), hence becoming stronger. As a result, wave-driven flows should emerge relatively easily in low-Prandtl-number fluids such as planetary cores made of liquid metal and stellar interiors [35,36], potentially affecting planetary and stellar dynamos [37] and magnetic reversals [38,39].

We now compare results of the full DNS model for the parameters of figure 1 (denoted by M 1 ) with results FIG. 2: Mean flow rms ūrms as a function of P r. The mean flow becomes stronger as P r decreases and is also more regular: the symbols' area is inversely proportional to the frequency bandwidth of ū, defined as the difference ∆f = f .9 -f .1 of the two frequencies f .1 and f .9 below and above which lies 10% of the mean-flow energy.

obtained from two reduced models (M 2 and M 3 ), described in figure 3A. The goal of the reduced models is to reproduce the evolution of the mean flow without resolving the convection. M 2,3 only solve the dynamics of the stable layer and are forced by prescribing values for the flow variables at its base (z N B ). If we force with exact time-varying values of (u , w ) and T from the full DNS, the evolution of ū is exactly reproduced in M 2 (not shown). Observations of real systems do not generally provide information about all variables at sufficient temporal and spatial resolution over long time periods, so we only use a subset of the full DNS data to force the reduced models. Specifically, here we expand the fluctuations u , w , T in M 2,3 in series of linear internal-wave modes, and we set their amplitudes such that the kinetic energy of each wave mode (defined by each wave's frequency ω and wave number k) matches the kinetic energy spectrum K(ω, k) obtained in the full DNS at z N B . We could have set the shape of the internal-wave spectrum by using theoretical predictions for the wave generation by turbulent convection [40,41], but that would preclude a comparison to M 1 , whose wave spectrum differs from, e.g., [40,41].

The reduced models only differ in how wave propagation away from the bottom boundary is solved. In M 2 , wave propagation is solved exactly by DNS of the Boussinesq equations, while in M 3 , a closed-form solution for the Reynolds stress is derived such that we only solve the 1D mean-flow (1). The analytical solution for the fluctuations (u , w ) in M 3 is obtained through WKB approximation, neglecting nonlinear wave-wave terms, mean-flow acceleration, and cross-interaction terms in the Reynolds stress (cf. details in Supplemental Material [START_REF]See Supplemental Material for details on the models, the effect of the Prandtl number, statistics and nonlinearities[END_REF]). We note that while M 2 is computationally cheaper than M 1 (resolution is 8 times smaller and time steps are ∼ 3 times larger), it remains significantly more demanding than M 3 , which is the only practical model for predicting the long-term dynamics of real systems (e.g. capturing the QBO in General Circulation Models). The goal of M 2 is to check approximations made in M 3 .

Figure 3B shows the temporal variations of the mean flow ū obtained in full DNS M 1 and in the two reduced models M 2 and M 3 . A large-scale oscillation is obtained in all three models, but the mean flow is stronger and the period is longer in the reduced models than in full DNS. Let us consider the characteristic amplitude of the mean-flow by its rms (ū rms ), and the characteristic period by taking the inverse of the peak frequency of its Fourier transform (T ū), which we average vertically between z = 0.8 and 1.1. We have: ūrms = 179 (M 1 ); 370 (M 2 ); and 415 (M 3 ); T ū = 0.125 (M 1 ); ≈ 0.33 (M 2 ); and ≈ 0.33 (M 3 ). Clearly, even when the wave propagation is solved exactly by DNS (M 2 ), the mean-flow dynamics is not reproduced quantitatively. In addition, the temporal variability of the mean flow obtained in full DNS is lacking in both reduced models.

The large discrepancy between the reduced models and DNS comes from the assumption that fluctuations on the lower boundary z = z N B of M 2,3 can be reconstructed from the time-averaged spectrum K using the linear wave relations for upward propagating plane waves. However, K can include contributions from overshooting plumes and some of the waves may be nonlinear. In our M 1 (M 2 ) simulations, the nonlinear terms have a typical magnitude of approximately 50% (10%) the linear terms just above z N B , suggesting overshooting convection in M 1 may be non-negligible at the interface. However, in the bulk of the stable region, this decreases to about 10% (5%), so the waves are in a weakly nonlinear regime (cf. details in Supplemental Material). Because the waves are weakly nonlinear, the energy transfer among waves does not affect the mean flow, explaining the agreement between M 2 and M 3 . Because forcing using the spectrum higher than z N B could in principle attenuate contributions from nonlinear convective motions, we have run additional simulations with different forcing heights (cf. Fig. S2 of the Supplemental Material [START_REF]See Supplemental Material for details on the models, the effect of the Prandtl number, statistics and nonlinearities[END_REF]): quantitative changes for the mean flow are obtained, but never lead to agreement with full DNS results. Importantly, the reconstruction of wave fluctuations from an energy spectrum neglects high-order statistics (higher than two), so statistics in the reduced models are Gaussian. However, intermittent events exist near the interface because the convection does not have a top-down symmetry and exhibits non-Gaussian statistics. In fact, the kurtosis of the fluctuations remains large, even in the wave field far from the interface z N B (see Fig. 3 of the Supplemental Material [START_REF]See Supplemental Material for details on the models, the effect of the Prandtl number, statistics and nonlinearities[END_REF]), suggesting that intermittency is a key component of wave generation. Intermittent intense wave events found in our DNS but neglected in the reduced models are typical of real systems. In the atmosphere, for instance, atmospheric waves sometimes propagate in the form of localized wave packets [42], such that wave intermittency can be non-negligible and has to be incorporated in reduced mean-flow models using stochastic processes [43].

In conclusion, the spontaneous generation and oscillation of a mean flow in our minimal, physical model is obtained for a wide range of parameters. In particular, we find that the mean flow becomes stronger as P r decreases (figure 2), which highlights the necessity to account for the real value of P r in stellar and planetary dynamical models. Evaluating the impact of wave-driven flows in natural systems is challenging. Indeed, we have shown here that reduced models do not yet predict the correct physics: tackling simultaneously the three-scale dynamics due to turbulence, waves, and mean flow, seems necessary. A major source of errors in reduced models comes from the approximations made in the types of waves excited by convection, even if the stably-stratified layer is forced with waves with the same kinetic energy spectrum as in full DNS. Our analysis suggests that imple-menting wave intermittency (through a boundary forcing scheme that would match the high-order moments of the DNS statistics), and disantengling non-wave contributions from the source spectrum are the next step forward and will be essential to improve the long-time predictive capabilities of low-order models.

The

Supplementary Information

SIMULATIONS DETAILS

Table I shows details (including typical time steps) of the simulations carried out for the full DNS model M 1 and reduced model M 2 with the open-source pseudo-spectral code Dedalus [START_REF] Burns | Dedalus: A Flexible Pseudo-Spectral Framework for Solving Partial Differential Equations[END_REF]. We use Chebyshev and Fourier modes in the z and x direction, respectively, and a 2-step implicit/explicit Runge-Kutta scheme for time integration. The CFL condition is 0.5. For reduced model M 3 , the mean-flow equation (with radiative damping from z = 1.35 to z = 1.5 as in M 1 & M 2 ) is solved using a second-order centered finite-difference scheme in z and an Adams-Bashforth/Crank-Nicolson scheme for time integration. The vertical resolution is δz = 0.00425 and time steps are δt = 0.0002.

To calculate the kinetic energy spectrum K, we record the values of the fluctuations u and w in the DNS for all x at a fixed interface depth z * and every δt = 5 × 10 -5 for one thermal time. In Fourier space, we thus have information about modes with wavenumbers ranging from π to 128π (discarding the 0th and Nyquist modes, L x = 2, n x = 256) and positive and negative angular frequencies ranging from 2π to 10 5 × 2π in absolute value (corresponding to retrograde and prograde waves, respectively), for a total of 2.56 × 10 6 modes. The highest-wavenumber modes put strong resolution constraints on the numerical simulation of the reduced models. Thus, we remove them from the dataset for simplicity and also because they carry so little energy that they do not affect the solution. Specifically, we remove all modes with wavenumber greater than 32π. Modes with angular frequencies greater than the buoyancy frequency N are not propagating modes so they are also removed from the dataset (in practice we remove all modes with angular frequency greater than 0.995N ). Finally, modes with angular frequencies lower than 0.005N are also removed because they carry little energy and because their treatment in M 3 requires special care since they reach critical layers easily. In a final step, we sum up every 8 successive energy bins in the frequency direction such that the final number of modes considered decreases to 32 × 788. Energy loss due to truncation of high-wavenumber, low-frequency and high-frequency modes when constructing the spectrum at height z * = 0.65, 0.68, 0.70 is small and reported in table II.

EFFECT OF DECREASING THE PRANDTL NUMBER

Figure 2 of the main text shows that the mean-flow rms becomes stronger and the flow becomes increasingly narrowbanded as the Prandtl number decreases. Because all physical effects in our model (waves, mean flow, convection) are fully coupled, changing P r changes not only the dissipation of the mean flow but also the convection and the wave generation and propagation. In order to provide a meaningful comparison of mean-flow emergence with changing P r, we adapted the other input parameters (Ra, S) of our DNS in such a way that the convection and the waves are not 2 in the main text). P r, Ra and S are input parameters chosen such that the x-averaged heat flux Q = -T z + wT , which is depth invariant, and x-averaged buoyancy frequency N (angular) are similar for all simulations. The domain size is fixed, i.e., L x × L z = 2 × 1.5, and the top temperature is also always T t = -43. Keeping N and T t constant fixes the interface depth to approximately the same value, as can be seen in supplementary figure 4. The CFL condition is set to 0.5 (δt gives the typical time step). For M 2 , P r and N are both input parameters. The intensification of the mean flow and narrowing of the mean-flow frequency spectrum with decreasing P r can be clearly seen from the full DNS results of ū shown as a function of t and z in supplementary figure 4. Notably, for P r ≥ 0.6, the QBO-like feature of the mean flow is weak or simply non-existent, while for P r = 0.1 the wave-driven mean flow penetrates deep inside the convective layer. The general trend shown in supplementary figure 4 suggests that for a given planetary or stellar state, defined by its heat loss and stable stratification strength, a QBO-like behavior will be more probable at low P r, which is the relevant limit for liquid iron and ionized plasmas.

EFFECT OF CHANGING THE BOTTOM BOUNDARY OF THE STABLE LAYER

Let us denote z * the bottom boundary of the stable layer in the reduced models M 2 and M 3 . In the main text, we show results for z * = z N B . When z * is varied (but not too different from z N B ), the mean flow exhibits similar patterns as in figure 3b M 2 -M 3 of the main text. However, this change does have a quantitative effect, as can be seen in supplementary figure 5 where we show the mean-flow rms ūrms and dominant period T ū as functions of z * . When z * increases, there is less energy available to force the mean flow (cf. table II) because internal waves are damped as they propagate upward. This leads to an increase of the mean-flow period as z * increases. The increase of T ū when energy of a broadband internal-wave spectrum decreases is in qualitative agreement with the results for monochromatic forcing [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF]. The effect of z * on ūrms is not monotonic, which is more surprising. In simulations with less energy (z * increases), one would expect that ūrms decreases. However, it is clear that ūrms and T ū remain much larger for all z * than what is obtained in full DNS (shown by the dashed lines). Thus, changing z * , while having a definite effect on the results, will not lead to values significantly closer to those in the full DNS. The differences between M 2 and M 3 are strongest for z = 0.65, with M 2 closer to the full DNS.

NON-GAUSSIAN STATISTICS IN FULL DNS

In the main text we suggest that wave intermittency may explain some of the discrepancies between the full DNS and the reduced models. In the reduced models, waves are forced from the bottom boundary assuming that they have the same amplitudes at all times. Thus, waves are not intermittent in the reduced models, and the statistics of the fluctuations are Gaussian (even in the bulk of the stable layer), which means that nonlinear effects do not make the flow statistics non-Gaussian. In the full DNS, however, the PDFs of the fluctuations are non-Gaussian both in the convection and the stably-stratified layers. Supplementary figure 6 shows the PDFs for w and u at z = z N B , and the standard deviation, skewness and kurtosis as functions of depth z. From the top row figure, it is clear that the tails of the PDFs have much larger values than what is obtained with a Gaussian distribution (shown by the dashed line) at z = z N B . From the bottom plots, we can see that the non-Gaussianity extends above the convective layer. While the standard deviation is much larger in the convective layer than in the stably-stratified layer (left plot), the positive skewness of vertical velocity (middle) remains relatively large in the entire domain. Positive skewness is expected in the convection since updrafts are more vigorous than downdrafts if there is an overlying stable layer (cf. [START_REF] Couston | [END_REF]). The kurtosis (right plot) is even more striking-the kurtosis of w is maximum close to z N B and decays only slowly toward the domain boundaries, where it reaches a kurtosis of 3, the result for Gaussian statistics. In summary, the large value of the kurtosis close to the interface best demonstrates that the process of wave generation and propagation away from z N B is characterized by intense, rare events (shown by the heavy tails in the top figure), and that wave intermittency should be included in parameterization schemes.

DETAILS OF THE REDUCED MODELS

In M 2 , we solve Eq. ( 2) of the main text via DNS with Dedalus for z ∈ [z * , L z ] assuming T = T + T , and using N 2 = -∂ z ρ = -P rRaS∂ z T = 4 × 10 8 , which corresponds to -P rRaS times the time and vertical average (from the interface height z * to L z ) of the temperature gradient in M 1 . Recall that overbar (prime) denotes x-average quantities (fluctuations).

For the reduced model M 3 , we solve Eq. ( 1) of the main text from z = z * to L z for the mean flow ū (which can be readily derived by taking the x average of Eq. ( 2) of the main text), with the fluctuations (u , w ) derived in closed form from the linear wave equation. Eq. (1) of the main text is solved numerically using finite differences and semi-explicit/implicit time stepping. At each timestep, (u , w ) are calculated from the mean flow profile, u(z), as described below; the Reynolds stresses are calculated assuming each wave only interacts with itself; then u(z) is updated by summing the effect of Reynolds stresses and viscosity.

The linearized wave equation can be written in terms of the stream function Ψ (u = -Ψ z and w = Ψ x ) under the mean-field approximation [45], as

∂ t + ū∂ x -∇ 2 ∂ t + ū∂ x -P r∇ 2 ∇ 2 -(ū zz + ūz ∂ z )∂ x Ψ = -N 2 ∂ 2 x Ψ. ( 4
)
As in M 2 , we take N 2 = -P rRaS Tz = 4 × 10 8 . Equation ( 4) can be solved analytically under WKB assumption, i.e. such that Ψ(x, z, t; ū) is obtained in closed form for all (x, z) and any mean-flow profile ū. Decompose the solution Ψ = Ψ + + Ψ -(± denotes waves going in the ±x direction) with Ψ ± of the form

Ψ ± = ψ ± e i(kx∓ωt) + * = A ± e z z * (iϕ±+ χ±)dz e i(kx∓ωt) + * (5)
with z * the interface height (taken as z N B in the main text), a small parameter, ω > 0, and * denotes the complex conjugate. Substituting (5) in (4), and assuming that dissipation terms are order O( ), t derivatives of ū are negligible, and z derivatives of ū are order O( ), yields

ϕ ± = ∓k 1 -2 ± 2 ± , χ ± = - 1 2 d dz ln |ϕ ± | - k 3 (1 + Pr) 2N 1 -2 ± 4 ±
, for upward-propagating waves, with ± = ω/N ∓ kū/N the normalized Doppler-shifted angular frequency. Convection generates a broad range of internal wave modes with different wavenumbers and frequencies such that the global streamfunction solution can be written as Ψ = k,ω ψ + e i(kx-ωt) + ψ -e i(kx+ωt) + * .

We decompose the Reynolds stress forcing of the mean flow as

-∂ z (w u ) = -∂ z (F + + F -) -∂ z C, (7)
where

F ± = - k,ω ikψ ± ∂ z ψ * ± + * = k,ω f ± (z) (8)
is the momentum fluxes due to the self interactions of prograde/retrograde waves and C is the momentum flux of all cross (i.e. non-self) interaction terms. Model M 3 neglects C and approximates the individual momentum fluxes using the WKB solution as

f ± = -kϕ ± (z = 0)|A ± | 2 exp   - Lz z * k 3 (1 + P r)dz N 4 ± 1 -2 ±   . (9)
Note that neglecting cross-interaction terms from the Reynolds stress is a customary approximation in GCMs, which makes the wave-driven forcing depends on time through ū only.

Internal waves have critical layers when their Doppler-shifted frequency goes to zero. In M 1 -M 2 , critical layers are regularized by dissipation or nonlinearity. In M 3 , when an internal wave approaches a critical layer, all of its remaining momentum is deposited following the work of [44]. We find that, in our configuration, setting the force due to internal waves at critical layers to zero yields almost the same results as depositing all of the momentum. This suggests that momentum deposition at critical layers is negligible compared to momentum deposition by viscous dissipation in the WKB model.

DETAILS OF THE FORCING AND BOUNDARY CONDITIONS IN THE REDUCED MODELS

The bottom forcing in the reduced models is constructed from the kinetic energy spectrum at the interface height z * extracted from full DNS, assuming that all motions correspond to linear internal waves (recall, again, that we take z * = z N B for the results of the main text). That is, we assume that the fluctuations can be written as a sum of wave modes of the form (sum runs over both prograde and retrograde modes)

(u , w , T ) = k,ω (U ± , W ± , T ± )e i(kx∓ωt+γ±z) + * (10)
where γ ± is the vertical wavenumber. We use the polarization relations of linear internal waves to express (U ± , W ± , T ± ) at z * in terms of K. K reads (hat denotes x, t Fourier transform)

K(ω, k) = 1 2 | û (z = z * )| 2 + | ŵ (z = z * )| 2 . ( 11
)
Using the linear dispersion relation γ ± = ∓k (N 2 -ω 2 )/ω 2 we obtain for the streamfunction amplitude (cf. ( 5)) and the other variables (ζ ± an arbitrary phase)

A ± = ωe iζ± kN K 2 , U ± = -iγ ± A ± , W ± = ikA ± , T ± = ±k Tz ω A ± . (12a)
Vertical boundary conditions are no-slip for the fluctuations at z = L z in M 2 , and ū = 0 at z = z * , L z in M 2 , M 3 . The no-slip bottom boundary condition for ū introduces errors, as ū = 0 at z = z * in full DNS. However, using other conditions does not lead to improvements of M 2 -M 3 results or introduces new tuning parameters. A stress-free boundary condition leads to unphysical jets forming close to the bottom boundary. A small damping layer (with, e.g., radiative damping coefficient D) with a free-slip boundary condition could potentially exert a sufficient drag on ū such that it reaches a physically-sound dynamical equilibrium, but it introduces the parameter D, which, with no prior study of its effect, could only be set in an ad-hoc manner.

QUANTIFICATION OF NONLINEAR EFFECTS IN THE WAVE REGION

The mean flow obtained in the different models may differ in part because of nonlinear effects in the wave region z ≥ 0.7. Therefore, in table III, we quantify the importance of nonlinear effects through the ratio N of the rms of the wave-wave to wave term |(u' ¯• ∇)u' ¯-(u' ¯• ∇)u' ¯| (overline denotes x average) divided by the rms of one of the dominant non-static force in the momentum equation, i.e. the dynamic pressure term |∇p |. The wave-mean to wave and wave-wave to mean nonlinear ratios N and N , given by the rms of |(ū ¯• ∇)u' ¯+ (u' ¯• ∇)ū ¯| and |(u' ¯• ∇)u' ¯| divided by the rms of |∇p |, are also reported in table III for completeness (other nonlinearities in the momentum equation are zero). Overall, N , N , and N decrease with height for z ≥ 0.7, and we provide values in table III at three different heights: close to the interface at z = 0.7, slightly above at z = 0.8, and in the bulk of the wave region at z = 1.0. For P r = 0.1, ū is strong, and we find N > N > N at z = 0.8, 1.0. For cases of P r ≥ 0.2, ū is weaker than for P r = 0.1 and the wave-wave to wave nonlinearity dominates; also, the relative importance of N compared to N increases as P r increases because the mean-flow weakens. For all M 1 simulations, N is relatively large (typically 50%) close to the interface, suggesting that wave-wave to wave nonlinear effects may be non-negligible; higher up, however, N is much smaller, typically 10% at z = 1.0, such that N is large at z = 0.7 mostly because there are traces of (strongly-nonlinear) convective motions at the base of the stable layer, and the waves above the interface are actually weakly nonlinear in all our simulations. We find that N is smaller in M 2 than in M 1 (P r = 0.2), even though the kinetic energy is the same at the interface. This finding is in agreement with the idea that nonlinear effects in M 1 are due to overshooting plumes at z = 0.7, and that waves have typically larger amplitudes (and more nonlinearity) in M 1 than in M 2 because they propagate in the form of intermittent intense wave packets in M 1 . N (N ) is larger (smaller) in M 2 than in M 1 , because ū is stronger but reverses more slowly in M 2 than in M 1 .

Model

FIG. 1 :

 1 FIG. 1: DNS results. (A) and (B) show snapshots of the vertical velocity field w at times t 1 = 2.21 and t 2 = 2.24, along with the mean flow ū (solid line) for z > z N B = 0.68 (with ū = 0 corresponding to x = 1). Vertical velocity patterns show convective motions in the lower part of the domain (z ≤ z N B ) and internal-wave motions in the upper part (z ≥ z N B). Note that energy propagates upward along wave crests, so crests toward the upper left (right) correspond to retrograde (prograde) waves. (C) shows the mean flow ū(t, z). The mean flow in the convective zone corresponds to the average of stochastic plumes emitted from the bottom boundary, hence reverses on a relatively rapid, convective time scale. In the stably-stratified layer, ū results from the nonlinear interaction of internal waves and oscillates on time scales ∼ 0.1, much longer than the buoyancy period ∼ π10 -4 . Simulation details and movies are available in Supplemental Material[START_REF]See Supplemental Material for details on the models, the effect of the Prandtl number, statistics and nonlinearities[END_REF].

FIG. 3 :

 3 FIG. 3: (A) Schematics of the DNS model M 1 and the two reduced models M 2 and M 3 . [1] We calculate the kinetic energy spectrum K of the fluctuations at height z N B obtained in M 1 . [2] The forcing (u , w , T ) is derived from K assuming that the fluctuations correspond to linear propagating internal waves. Propagation of the waves is solved [3] via DNS of the Boussinesq equations in M 2 , but is derived analytically [4] in M 3 under WKB approximation. Thus, in M 3 , [5] we only need to solve for the mean-flow equation. As in full DNS, we use a damping layer for 1.35 < z < 1.5, and boundary conditions for the mean flow are no slip. (B) ū over one thermal time obtained for each model shown in (A). Physical parameters are as in figure 1. Note that the colormap has been changed in figure 3B compared to figure 1C to highlight differences between M 1 and M 2,3 .

  full K(z = 0.65) truncated K(z = 0.65) full K(z = 0.68) truncated K(z = 0.68) full K(z = 0.70) truncated K(z = 0

FIG. 4 :

 4 FIG.4: Mean flow ū obtained in DNS for different Prandtl numbers, shown as a function of z over one thermal time after reaching statistically steady state. Time on the x axis is the reduced time with t * = 0.6, 2.0, 4.0, 2.0, 2.0, 2.0 for P r = 0.1, 0.2, 0.3, 0.6, 1.0, 3.0. We use the same color axis for ū in all plots in order to highlight the differences of the mean-flow amplitude in the stable region for all P r.

FIG. 5 :

 5 FIG. 5: Variations of the dominant mean-flow period (left plot) and mean-flow rms (right) averaged from z = 0.8 to z = 1.1, with the height z * , the bottom of the stable layer. Results are shown as open circles for M 2 and as filled upper triangles for M 3 . The full DNS results for ūrms and T ū are shown by the dashed lines.

FIG. 6 :

 6 FIG. 6: (top) The normalized PDFs of w and u at z = z N B = 0.68 with solid lines, and the Gaussian distribution with the dashed line. The means of w and u are 0. The bottom row shows the standard deviation (left), skewness (middle) and kurtosis (right) of w and u recorded over one thermal time for the DNS results presented in the main text at 17 different depths (shown by symbols).

TABLE I :

 I Parameters of DNS with varying Prandtl number P r (figure

TABLE II :

 II Total kinetic energy spectrum K for the full DNS dataset and after modal truncation (see the text) at heights z = 0.65, 0.68, 0.70.

strongly affected. Specifically, we have set Ra and S such that the x-averaged vertical heat flux Q is constant, which fixes the energy available to convection, and such that the buoyancy frequency N is constant, so wave propagation is not affected. All physical parameters are listed in table I.

TABLE III :

 III Quantification of nonlinear effects in M 1 and M 2 simulations. All numbers are in %. See the text for the definition of N , N , N .