2. Phenomenology

To gain insight into the problem, it is relevant to consider the prototypical case where the stirring velocity field has a constant kinetic energy $\mathcal{E}_{k}$ (per unit mass). In that case, the long-time limit of the cumulative mixing efficiency $\eta$ only depends on $\mathcal{E}_{k}$ and the initial density field configuration. For instance, the background potential energy (per unit volume, scaled by a reference density $\rho _{0}$ ) for an initially linearly stratified density field $\rho (z, t=0)=-N^{2}\rho _{0}z/g$ (with $z$ the vertical coordinate, $N$ the buoyancy frequency and $g$ the acceleration of gravity) and a fully mixed final state (characterised by a single density level $\langle \rho \rangle := \int _{0}^{H}\rho (z,t=0)\textrm{d}z / H = -N^{2}\rho _{0}H/(2g)$ with $H$ the depth of the water column) is

(2.1) \begin{equation} \mathcal{E}_{b}(t=0)=\frac {g}{\rho _{0}H}\int _{0}^{H}\rho (z,t=0)z\,\textrm{d}z=-\frac {1}{3}N^{2}H^{2}, \end{equation}

and

(2.2) \begin{equation} \mathcal{E}_{b}(t \to +\infty )=\frac {g}{\rho _{0}H}\int _{0}^{H}\langle \rho \rangle z\,\textrm{d}z=-\frac {1}{4}N^{2}H^{2}, \end{equation}

respectively. Hence,

(2.3) \begin{equation} \eta (t \to +\infty ) = \frac {1/12}{1/12 + 1/Ri} \quad \text{with} \quad Ri=\frac {N^{2}H^{2}}{\mathcal{E}_{k}}. \end{equation}

In the two-layer configuration depicted in figure 1, a similar computation gives

(2.4) \begin{equation} \eta (t \to +\infty ) = \frac {1/8}{1/8 + 1/Ri} \quad \text{with} \quad Ri= \frac {gH}{\mathcal{E}_{k}}\times \frac {\rho _{0} - \rho _{1}}{\rho _{0}}. \end{equation}

The above-mentioned long-time limits of the mixing efficiency $\eta$ are decreasing functions of the $\mathcal{E}_{k}$ : the more energy needing to be input into the system to fully homogeneise it, the lower the mixing efficiency. Note that for $Ri \lt 1$ , the kinetic energy in the flow is sufficient to overturn the full water column; for $Ri \gt 1$ , it is not. Equation (2.4) justified the constant kinetic energy scenario; the long-time limit of the mixing efficiency is fixed and we here focus on its time evolution towards that limit.

Two processes control the time evolution of $\eta$ towards its long-time limit: stirring that reorganises the density field without changing $\eta$ , and mixing that smooths the density field, increases the background potential energy $\mathcal{E}_{b}$ and hence $\eta$ . We study these two processes in isolation and then in combination.

Figure 1. Evolution of a two-layer density field in the absence (top) or presence (bottom) of mixing. The density level $\rho _{0}$ and $\rho _{1}\lt \rho _{0}$ are represented in green and white, respectively; shades of green correspond to intermediate density levels. As turbulent fluctuations act against stratification to stir the box and vertically disperse fluid parcels with different densities, they create a steady-state density field characterised by a well-scrambled region where dense and light fluid coexist ((a) $\rightarrow$ (b)). This steady state is reached within a time $t_{R}$ . Thick black lines in panel (b) correspond to the steady-state horizontally averaged density field derived by Venaille, Gostiaux & Sommeria (Reference Venaille, Gostiaux and Sommeria2017) ( $\overline {\rho }\propto \tanh {(Riz/H)}$ ) for various values of the Richardson number $Ri$ . Note that in the two-layer configuration, the well-scrambled region has a thickness $L_{S}/H \sim 1/Ri$ (rather than $L_{S}/H \sim 1/\sqrt {Ri}$ as in the linearly stratified case). The blue curves correspond to the probability of finding a given density level at a given point in space $p(\rho , z)$ . The effect of mixing is to smooth density differences between neighbouring fluid parcels ((c) $\rightarrow$ (d)) within a time $t_{M}$ . Here, we model the effect of mixing as an averaging process in density space, i.e. a (weighted) convolution in probability space. In this work, we will assume that the volumes of the fluid parcels mixing into each other are equal ( $V_{0}=V_{1}$ ) so that $\alpha :=V_{0}/(V_{0}+V_{1})=1/2$ .

In the absence of mixing ( $\kappa =0$ ), turbulent fluctuations (stirring) and buoyancy forces (restratification) compete so that the vertical displacement of fluid parcels with a given density content (that is constant since $\kappa = 0$ for now) is constrained (Kimura & Herring Reference Kimura and Herring1996). In the initially linearly stratified configuration, the furthest (dimensionless) vertical distance such fluid parcels can travel away from their initial position is $L_{S}/H \sim \sqrt {\mathcal{E}_{k}/(N^{2}H^2)} = \sqrt {1/Ri}$ (a result that arises from a force balance between inertia and buoyancy on fluid parcels; see Petropoulos, de Bruyn Kops & Caulfield (Reference Petropoulos, de Bruyn Kops and Caulfield2025)). In the two-layer configuration of figure 1, we have $L_{S}/H \sim 1/Ri$ (Venaille et al. Reference Venaille, Gostiaux and Sommeria2017). The weaker the stratification, or the more energy there is in the system, the further away fluid parcels can vertically travel. This has important implications for mixing, as the ‘screening length’ $L_{S}$ inevitably constrains the density levels that fluid parcels encounter along their stochastic path and hence mix into; for $L_{S}/H \ll 1$ (i.e. $Ri \gg 1$ ), the density levels seen by fluid parcels are close to the density they carry and we expect the route to homogeniety to be slow. By stirring the density field with constant kinetic energy, we fix this parameter.

As a consequence of the competition between stirring and restratification, in the diffusiveless limit $\kappa = 0$ , the probability $p(\rho , z)$ of finding a density level $\rho$ at depth $z$ reaches (maximum disorder) a steady-state limit (Venaille et al. Reference Venaille, Gostiaux and Sommeria2017) characterised by ‘well-scrambled’ (or ‘overturned’) regions of dimensionless depth $L_{S}/H$ (figure 1). We hypothesise here that the characteristic time scale to reach this equilibrium state (i.e. to overturn layers of depth $L_{S}$ ) is the time for a fluid parcel to come back to its equilibrium position due to restoring buoyancy forces after being vertically disturbed. Because velocity gradients typically stretch fluid parcels into thin and elongated structures, or lamellae, with characteristic thickness below the scale at which viscosity starts to matter (Duplat et al. Reference Duplat, Innocenti and Villermaux2010; Villermaux Reference Villermaux2019), this time is not the inverse of the buoyancy frequency, but rather the resetting time $t_{R} \sim Re_{b}/\gamma \sim \sqrt {Re_{b}}/N$ with $\gamma$ the average stretching rate, and $Re_{b} := \epsilon /(\nu N^{2})$ the buoyancy Reynolds number and $\epsilon$ the dissipation rate of turbulent kinetic energy (Petropoulos et al. Reference Petropoulos, de Bruyn Kops and Caulfield2025). Here, we used the estimate $\gamma \sim \sqrt {\epsilon /\nu }$ for the (average) stretching rate (Batchelor Reference Batchelor1959). The viscous correction is crucial. It takes into account that the settling dynamics of the lamellae towards their neutrally buoyant positions is dictated by a balance between viscous drag and buoyancy forces due to the small characteristic scales involved (Petropoulos et al. Reference Petropoulos, Caulfield, Meunier and Villermaux2023, Reference Petropoulos, de Bruyn Kops and Caulfield2025). Note that in the two-layer configuration of figure 1, this balance gives $t_{R} \sim \gamma \mathcal{E}_{k}\rho _{0}^{2}/[g (\rho _{0} - \rho _{1})]^{2} \sim \sqrt {\epsilon /\nu } \times \mathcal{E}_{k}\rho _{0}^{2}/[g (\rho _{0} - \rho _{1})]^{2}$ .

The screening length scale $L_{S}$ and resetting time $t_{R}$ can be used to define a ‘dispersivity’ $D := L_S^{2}/t_{R}$ characterising the dispersive effect of turbulent fluctuations in the presence of stratification (in the diffusiveless limit $\kappa = 0$ ).

We now consider the effect of mixing ( $\kappa \neq 0$ ). As turbulent fluctuations disorder the density field, small-scale density structures (the lamellae) with different density content will encounter each other and mix; this happens on the mixing time scale $t_{M} \sim \gamma ^{-1}\mathcal{F}(Pr) \sim \sqrt {\nu /\epsilon }\mathcal{F}(Pr)$ (Petropoulos et al. Reference Petropoulos, Caulfield, Meunier and Villermaux2023, Reference Petropoulos, de Bruyn Kops and Caulfield2025). Then, if $t_{M} \lt t_{R}$ , we expect $\eta$ to reach its long-time limit relatively fast; as soon as parcels are vertically disturbed, they mix with their surroundings. Conversely, if $t_{M} \gt t_{R}$ , turbulence first stirs the density field towards its maximum disorder equilibrium state before mixing and the time evolution of $\eta$ is slower; in particular, mixing is not efficient in the sense that vertically disturbed fluid parcels do not have time to mix before returning to their neutrally buoyant position. More precisely, we hypothesise the existence of the following regimes in the ( $Ri$ , $t_{M}/t_{R}$ ) parameter space (summarised in figure 2).

1. (i) Full overturn of the water column, i.e. $Ri\lt 1$ . In that case, there is enough energy to overturn the full water column ( $L_{S} \gt H$ ) and, hence, for fluid parcels to potentially encounter all the density levels present in the column along their stochastic paths. If the mixing time is larger than the time for turbulent fluctuations to overturn the entire water column (i.e. $t_{M} \gt H^{2}/D$ or, equivalently, $Ri \lt t_{M}/t_{R}$ ), mixing is mixing-limited. We expect the water column to be fully mixed after a time of order $t_{R}+t_{M}$ , i.e. when scaled by $t_{R}$ , of order $t_{M}/t_{R}$ in the limit $t_{M} \gg t_{R}$ . The turbulent fluctuations first overturn the water column, then mixing happens at once. In contrast, if $t_{M} \lt H^{2}/D$ , mixing of the water column is limited by dispersion. Mixing happens before turbulent fluctuations have time to overturn the water column. We expect that, after a mixing time $t_{M}$ , well-mixed layers of depth $\sqrt {Dt_{M}}$ are formed (this corresponds to the vertical excursion of fluid parcels due to the dispersivity $D$ within a mixing time). There are $H/\sqrt {Dt_{M}}$ such layers to mix with each other. Hence, the time to fully mix the water column is expected to be of order $t_{M} \times H/\sqrt {Dt_{M}}$ , i.e. $\sqrt {Ri\times t_{M}/t_{R}}$ when scaled by $t_{R}$ .

2. (ii) Partial overturn of the water column, i.e. $Ri\gt 1$ . In that case, the energy in the system can only overturn layers of depth $L_{S}\lt H$ . In the case where mixing happens before turbulent fluctuations can overturn the depth $L_{S}$ of the water column that is accessible energetically (i.e. $t_{M} \lt L_{S}^{2} / D$ or, equivalently, $t_{M}/t_{R} \lt 1$ ), mixing is dispersion limited. The above-mentioned argument holds and we expect the water column to be fully mixed after a (dimensionless) time of order $\sqrt {Ri \times t_{M}/t_{R}}$ . In contrast, if $t_{M}/t_{R} \gt 1$ , we expect that after $t_{R}$ , there are $H/L_{S}$ ‘well-scrambled’ layers (in that regime, mixing has no time to act before fluctuations overturn layers of depth $L_{S}$ ) that need to be mixed internally and with each other. Hence, the time to mix the entire water column is expected to be of order $t_{M} \times H/L_{S}$ , or when scaled by $t_{R}$ , $\sqrt {Ri} \times t_{M}/t_{R}$ .

Figure 2. (a) Regime diagram in $(Ri, t_{M}/t_{R})$ space. (b) Time evolution of the dissipation rate of turbulent kinetic energy $\epsilon$ . (c) Time evolution of the ratio of time scales $t_{M}/t_{R}$ . The (dimensionless) time scale governing how fast the water column is fully mixed is presented as text insets. Symbols correspond to simulation data: $Ri_L=0.1$ , $H=1$ ; $Ri_{L}=0.5$ , $H=1$ ; $Ri_{L}=0.1$ , $H=\sqrt {10}$ ; $Ri_L=1$ , $H=1$ ; $Ri_L=0.5$ , $H=\sqrt {10}$ ; $Ri_L=1$ , $H=\sqrt {10}$ ; $Ri_{L}=10$ , $H=10$ ; $Ri_L=20$ , $H=1$ ; $Ri_L=30$ , $H=1$ and $Ri_L=50$ , $H=1$ .

The different elements described previously can be summarised into the following equation for the evolution of the probability $p(\rho ,z)$ of finding a density level $\rho$ at depth $z$ :

(2.5) \begin{equation} \partial _{t}p = \underbrace {D\partial _{z}^{2}p}_{\text{dispersion}} + \underbrace {\frac {Dg}{\beta \rho _{0}}\partial _{z}[(\rho -\overline {\rho })p]}_{\text{restratification}} + \frac {1}{t_{M}}\underbrace {(-p + p_{\alpha } \ast p_{1-\alpha })}_{\text{mixing}}. \end{equation}

Here, $\overline {\rho }(z):=\int \rho p(\rho ,z)\,\textrm{d}\rho$ is the instantaneous mean density profile and $p_{\alpha }(\rho , z) = p(\rho /\alpha , z)/\alpha$ . The convolution operation $\ast$ is in density space.

The dispersion and restratification parts of (2.5) follow from work by Venaille et al. (Reference Venaille, Gostiaux and Sommeria2017). The constant $\beta$ is chosen so that, in the diffusiveless limit $t_{M}/t_{R} \to +\infty$ , (2.6) has a steady-state solution that exhibits a characteristic length scale $L_{S}$ , i.e. $D/L_{S}^{2} \sim Dg/(\beta L_{S})$ . The dispersion operation models turbulent fluctuations that broaden the distribution of density levels found at a given height $z$ . Restratification towards the instantaneous mean density profile $\overline {\rho }$ models the effect of buoyancy that increases the probability of finding a fluid parcel close to its neutrally buoyant (with respect to $\overline {\rho }$ ) position. Mixing is modelled as a convolution operation, translating the fact that mixing averages the density of neighbouring fluid parcels within a mixing time (figure 1). Importantly, by construction, this operation preserves the average density and its stationary solution is the Dirac measure $\delta (\rho - \langle \rho \rangle )$ (see Appendix A). Here, we assume that the weight $\alpha = 1/2$ , i.e. that the different fluid structures mixing into each other have approximately the same volume.

When scaling time $t$ by the resetting time $t_{R}$ , height $z$ by the depth of the water column $H$ and the density $\rho$ by the reference density $\rho _{0}$ , (2.5) reads

(2.6) \begin{equation} \partial _{t}p = Ri^{-s}\partial _{z}^{2}p + Ri^{-s/2}\partial _{z}[(\rho -\overline {\rho })p] + \frac {t_{R}}{t_{M}}(-p + p_{\alpha } \ast p_{1-\alpha }). \end{equation}

In the linearly stratified configuration, $s=1$ and $t_{M}/t_{R} \sim \nu N^{2}\epsilon ^{-1}\mathcal{F}(Pr)$ . In the two-layer system depicted in figure 1, $s=2$ and $t_{M}/t_{R} \sim \nu g^{2}(\rho _{0}-\rho _{1})^{2}/(\epsilon \mathcal{E}_{k}\rho _{0}^{2})\mathcal{F}(Pr)$ . Note that mass conservation is ensured by the no-flux boundary conditions $Ri^{-1}\partial _{z}p + Ri^{-1/2}(\rho -\overline {\rho })p = 0$ for all $\rho$ at $z=0,1$ (see Appendix A). This boundary condition and (2.6) model the time evolution of $p$ through stirring and mixing. The background potential energy and, hence, the mixing efficiency $\eta$ can then be estimated from $p$ using the method of Tseng & Ferziger (Reference Tseng and Ferziger2001).

3. Comparison with direct numerical simulations

We now compare the theoretical scalings described previously with three-dimensional direct numerical simulations. In this section, dimensional quantities are denoted by a hat $\hat {\boldsymbol{\cdot }}$ . Dimensionless quantities are left unchanged. We use the spectral element solver Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020) to solve the (dimensionless) incompressible Navier–Stokes equation under the Boussinesq approximation

(3.1) \begin{equation} \begin{cases} \partial _{t}\boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} = -\boldsymbol{\nabla }p + Re_{L}^{-1}{\nabla} ^{2}\boldsymbol{u} - Ri_{L}\rho \boldsymbol{e}_{z} + A\boldsymbol{u}, \quad \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} = 0, \\ \partial _{t}\rho + \boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{\nabla }\rho = (Pr Re_{L})^{-1}{\nabla} ^{2}\rho , \end{cases} \end{equation}

with stress-free boundary conditions and no density flux at $z=0, H$ , and periodic boundary conditions in the horizontal. The forcing term $A\boldsymbol{u}$ is tuned at each time-step so as to have a bulk kinetic energy $\mathcal{E}_{k} = 1$ at all times, following the method of Yi & Koseff (Reference Yi and Koseff2022). In particular, $A$ satisfies

(3.2) \begin{equation} A = \frac {\epsilon + \mathcal{B} - \tau ^{-1}(\mathcal{E}_{k} - 1)}{2\mathcal{E}_{k}}, \end{equation}

with $\mathcal{B}$ the buoyancy flux. This ensures that a constant kinetic energy $\mathcal{E}_{k}=1$ is reached within $\tau$ units of time. Here, $\tau =0.1$ (in dimensionless, simulation time units).

We initialise the simulation with a linearly stratified density field with dimensional buoyancy frequency $\hat {N}$ and a divergent-free, isotropic velocity with dimensional energy spectrum (Passot & Pouquet Reference Passot and Pouquet1987)

(3.3) \begin{equation} \hat {E}(\hat {k}) = \frac {32 \hat {E}_{0}}{\hat {k}_{0}}\sqrt {\frac {2}{\pi }}\left (\frac {\hat {k}}{\hat {k}_{0}}\right )^{4}\exp \left [-2\left (\frac {\hat {k}}{\hat {k}_{0}}\right )^{2}\right ]. \end{equation}

The initial kinetic energy $\hat {E}_{0}$ and wavenumber with maximum energy $\hat {k}_{0}$ determine a velocity and length scale $\hat {U}:=\sqrt {\hat {E}_{0}}$ and $\hat {L}:=1/\hat {k}_{0}$ , respectively. We use the dynamic pressure $\hat {\rho }_{0}\hat {U}^{2}$ to scale pressure and $\hat {L}\hat {\rho }_{0}\hat {N}^{2}/\hat {g}$ to scale density (with $\hat {\rho }_{0}$ a reference density). Time is scaled advectively so that $Re_{L} := \hat {U}\hat {L}/\hat {\nu }$ and $Ri_{L} := \hat {N}^{2}\hat {L}^{2}/\hat {U}^{2}$ . The dimensionless horizontal box size is $10 \times 10$ . Simulations are run until the density field is homogeneous. We fix $Re_{L} = 100$ , and choose the value of $Ri_{L}$ and $H:=\hat {H}/\hat {L}$ to target specific values of $Ri = \hat {N}^{2}\hat {H}^{2}/\hat {\mathcal{E}_{k}} =Ri_{L}H^{2}/\mathcal{E}_{k}$ (with $\hat {\mathcal{E}_{k}}$ scaled by $\hat {U}^{2}$ ). The value of $\hat {\epsilon }$ and, hence, $t_{M}/t_{R} \sim \hat {\nu } \hat {N}^{2}\hat {\epsilon }^{-1}\mathcal{F}(Pr) = Ri_{L}\langle \partial _{i}u_{j}\partial _{i}u_{j}\rangle _{V}^{-1}\mathcal{F}(Pr)$ (with $\langle \boldsymbol{\cdot }\rangle _{V}$ being a volume average; here, $\hat {\epsilon }$ is scaled by $\hat {\nu } \hat {U}^{2}/\hat {L}^{2}$ ) arises as a property of the flow, and is in fact time dependent (see figures 2 b and 2 c) since no statistical equilibrium can be reached as no background stratification is imposed. More precisely, as velocity gradients develop in the flow, $\epsilon$ grows. When the density field is fully mixed and $\mathcal{B}=0$ , the value of $\epsilon$ is essentially fixed by the forcing parameter $A$ (see (3.2)). The forcing chooses to destroy small-scale velocity fluctuations and drives the water column at just the right speed to have $\mathcal{E}_{k}=1$ . As a result, $\epsilon$ decreases (see figure 2 b). We will assume, admittedly somewhat cavalierly, that the parameter $t_{M}/t_{R}$ is well characterised by its time average. We fix $Pr=1$ so that $\mathcal{F}(Pr) \simeq 1$ . The simulation parameters are summarised in figure 2.

Figure 3. (a) Time evolution of the cumulative mixing efficiency $\eta$ (symbols) with prediction from model (2.6) (solid lines). In panels (b) and (c), time has been rescaled using predictions of § 2in the ‘fast’ mixing () and ‘slow mixing’ () regimes. (d) Same as panel (a) for the time-evolution of the density variance, scaled by its initial value.

The time evolutions of $\eta$ are presented in figure 3 with the corresponding prediction from the model (2.6) (numerically solved using 20 depth levels and 100 density levels). Note that the simulation time is multiplied by $Ri_{L}\langle \partial _{i}u_{j}\partial _{i}u_{j} \rangle _{V}^{-1/2}$ to be consistent with the dimensionless time used in (2.6). We found that the best fits occur for $t_{M}/t_{R} \simeq 5.5\hat {\nu } \hat {N}^{2} \hat {\epsilon }^{-1} = 5.5/Re_{b}$ . When time is rescaled by the algebraic combination of $Ri$ and $t_{M}/t_{R}$ suggested in § 2, the data collapse on a single curve, suggesting that these times indeed control the evolution of the density field from its initial to a fully mixed state.

Another important measure of mixing is the time evolution of the density histogram variance, characterising the diversity of density levels still present in the flow as it mixes. In figure 3(d), we present the time evolution of the density variance as it decreases from its initial value (characterising the spread of density levels present in the initial configuration of the density field) to $0$ (i.e. fully homogeneous density field, no variance). Again, the model (2.6) appears to capture the simulation trends, suggesting that the modelling assumptions used in this work are relevant to describe mixing.