2.1. Equations of motion

Throughout this study, we designate the two density-changing components to be ‘salinity’ and ‘temperature’ for respectively the least and the most diffusing components. We consider a layer of fluid with thickness $d$ enclosed between two parallel infinite walls that are maintained at constant temperature and salinity. The temperature and salinity differences between the top and the bottom walls are assumed to be ${\rm \Delta} T$ and ${\rm \Delta} S$, respectively. The control parameters of the problem are the Prandtl, Lewis, thermal and haline Rayleigh numbers, defined as

(2.1a–d)\begin{equation} Pr = \frac{\nu}{\kappa_T},\quad Le = \frac{\kappa_T}{\kappa_S},\quad R_T = \frac{g\alpha{\rm \Delta} T d^3}{\nu \kappa_T},\quad R_S = \frac{g\beta{\rm \Delta} S d^3}{\nu \kappa_S}, \end{equation}

where $\nu$ is the kinematic viscosity, $\kappa _T$ is the thermal diffusivity, $g$ is the acceleration due to gravity, and $\alpha$ and $\beta$ are the thermal and haline linear expansion coefficients. From these, the density ratio $R_{\rho }$ is defined as

(2.2)\begin{equation} R_{\rho}=Le\frac{R_T}{R_S}. \end{equation}

The fluid motion is governed by the Boussinesq equations that take the following dimensionless form (Paparella & von Hardenberg Reference Paparella and von Hardenberg2012):

(2.3a)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}

(2.3b)\begin{gather}\frac{\partial\boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u} ={-}\boldsymbol{\nabla} p +Pr Le\left(\nabla^2\boldsymbol{u}+R_SB \hat{\boldsymbol{z}}\right), \end{gather}

(2.3c)\begin{gather}\frac{\partial T}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})T = Le\nabla^2 T, \end{gather}

(2.3d)\begin{gather}\frac{\partial S}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})S = \nabla^2 S. \end{gather}

In this paper, we consider the two-dimensional case, thus $\boldsymbol {u}=(u, w)$ is the velocity field of the fluid that occupies the dimensionless domain $\boldsymbol {x}=(x,z)\in [0, L_x]\times [0, 1]$; $p$ is the pressure; $T$ and $S$ are the temperature and salinity fields, whose value ranges from 0 at the bottom wall to 1 at the top wall; $B= R_{\rho }T-S$ is the buoyancy field and $\hat {\boldsymbol {z}}$ is the upward unit vector. The vorticity field associated with $\boldsymbol {u}$ is defined to be $\omega :=\partial _{z}u-\partial _{x}w$. All variable fields are assumed to be laterally periodic. Free-slip conditions are enforced on momentum at both walls. The stationary base state of (2.3a)–(2.3d) reads

(2.4)\begin{equation} \left(\boldsymbol{u}_b, T_b, S_b\right) = \left((0, 0), z, z\right). \end{equation}

The solutions of (2.3a)–(2.3d) develop thermal and haline boundary layers near the top and bottom walls. In the interior, depending on the value of the parameters in (2.1a–d), the horizontally averaged temperature and salinity profiles may either form a slowly evolving staircase-like profile or reach a statistically steady state characterized by constant vertical gradients $G_T$ and $G_S$, with $G_S < G_T < 1$. In this case, the strength of the convection is best quantified by the local density ratio $R_{\rho \ {loc}}$, defined as

(2.5)\begin{equation} R_{\rho\ {loc}} = \frac{G_T}{G_S}R_{\rho}. \end{equation}

We define the anomaly of a scalar field $({\cdot })$ with respect to its horizontal average as

(2.6)\begin{equation} ({\cdot})':= ({\cdot})-\overline{({\cdot})}:=({\cdot})- \frac{1}{L_x} \int_{0}^{L_x} ({\cdot})\,\mathrm{d}\kern0.7pt x, \end{equation}

in which $\overline {({\cdot })}$ is the horizontal average of $({\cdot })$. The buoyancy anomaly is thus

(2.7)\begin{equation} B^\prime = R_{\rho}T^\prime-S^\prime. \end{equation}

Because $\bar {w}=0$, only the buoyancy anomaly flux $B^\prime w$ contributes to the total advective buoyancy flux across any horizontal section of the domain.

We define a Lagrangian particle as a geometric point whose velocity coincides with that of the fluid. That is, a Lagrangian particle sweeps a curve $\boldsymbol {x}_0\mapsto \boldsymbol {x}(\boldsymbol {x}_0; t)$ uniquely defined by

(2.8a,b)\begin{equation} \frac{{\rm d}\kern0.7pt \boldsymbol{x}(t)}{{\rm d} t} = \boldsymbol{u}(\boldsymbol{x}(t), t),\quad \boldsymbol{x}(\boldsymbol{x}_0; 0) = \boldsymbol{x}_0, \end{equation}

where $\boldsymbol {u}(\boldsymbol {x}, t)$ is the fluid velocity at the current particle position (as given by (2.3a)–(2.3d)). We designate a Lagrangian particle's buoyancy, salinity and temperature anomalies to be the same as those of the fluid's at the particle's current location.

2.2. Numerical set-up

Equations (2.3a)–(2.3d), recast in the vorticity-streamfunction formulation, were solved numerically with a Galerkin pseudo-spectral code. The numerical method uses the Fourier basis in the horizontal direction and the sine basis in the vertical. Time evolution is implemented through a third-order Adams–Bashforth scheme. The initial condition consists of the conductive solution (2.4) with a small white noise perturbation added to the salinity field. We monitor the temperature and salinity fluxes computed in the bulk of the domain (defined as the middle third of it) and at the boundary of the domain. During the initial transient, the fluxes in the bulk are larger (in absolute value) than the fluxes at the boundary. When the transient is over and the boundary layers become fully established, the fluxes in the bulk and at the boundary match each other, and remain nearly constant in time. To check that no thermohaline staircases are formed, we extend the Eulerian part of the simulations for at least five times the duration of the initial transient (see supplementary figure S7 available at https://doi.org/10.1017/jfm.2023.311). During this time, we check that a statistically steady state is maintained where the horizontally averaged temperature, salinity and buoyancy, far from the boundaries, keep a constant vertical gradient. We then seed at uniformly random positions within the entire computational domain a number of Lagrangian particles equal to the number of Eulerian grid points. The time step (${\rm \Delta} t$) used in the simulations was chosen to maintain a Courant–Friedrichs–Lewy (CFL) number (based on the vertical velocity) of approximately $0.2$. The Lagrangian data are saved every 100 time steps. Each Lagrangian simulation lasts for a time interval in which a particle travelling at the CFL speed would cover five times the height of the computational domain. The motion of the Lagrangian particles is tracked with the same time integrator used for the Eulerian part. The velocity of the particles is computed by taking the derivative of a bicubic interpolator of the streamfunction. This method ensures the divergenceless of the interpolated velocity field. The parameters of the simulations used in this study are presented in table 1. We consider two groups of simulations: one where the density ratio is fixed to $R_\rho =1.2$ and the haline Rayleigh number is $R_S=10^{12},\ 3.16\times 10^{12},\ 10^{13},\ 3.16\times 10^{13},\ 10^{14}$, and one in which the haline Rayleigh number is fixed to $R_S=10^{12}$ and the density ratio is $R_\rho =2.8, 2.4, 2.0, 1.6, 1.2$. In all simulations, we fixed $Pr =10$ and $Le = 3$. The simulations with fixed density ratio explore a regime very far from marginality, where fingering convection is quite vigorous and is characterized by large fluxes and fast-moving blobs, but not vigorous enough to lead to staircase formation. The simulations with fixed haline Rayleigh number are intended to investigate whether the blob-dominated regime breaks down when the density ratio approaches the marginally stable value $R_\rho =Le$.

Table 1. Parameters of simulation runs: haline ($R_S$) and thermal ($R_T$) Rayleigh numbers, domain width ($L_x$), mesh size ($N_x\times N_z$), time step (${\rm \Delta} t$) and local density ratio $R_{\rho \ {loc}}$ at equilibrium. We fixed $Pr =10$ and $Le = 3$ in all simulation runs. Note that the simulations with $R_S=10^{12}$ have a density ratio $R_\rho =2.8, 2.4, 2.0, 1.6, 1.2$, from top to bottom; the other simulations all have $R_\rho =1.2$.

2.3. Lagrangian legs and leg statistics

We partitioned the trajectory of every Lagrangian particle into a set of shorter fragments that we refer to as the legs. We define a leg to be a maximal section of a trajectory with strictly increasing or decreasing elevation, i.e. the $z$-coordinate along a leg is monotone. The monotonicity alternates between adjacent legs. For every leg, $l$, we define the signed leg length, $\zeta$ of $l$, to be the difference between $z$-coordinates of its terminal and initial points. Thus, legs with a positive length correspond to up-going particles, and those with a negative length correspond to down-going particles. In the supplementary figure S1, we present a visual illustration of a particle trajectory and of the legs of which it is made.

We define a leg's properties (buoyancy anomaly, vertical velocity, flux of buoyancy anomaly, etc.) to be the mean of those quantities along the leg. Specifically, if a leg begins at time $t_s$, ends at time $t_e$ and describes the curve $\boldsymbol {x}$, then we define its property $g_l$ as

(2.9)\begin{equation} g_l := \frac{1}{t_e - t_s} \int_{t_s}^{t_e} g(\boldsymbol{x}(t), t)\,{\rm d} t\quad \mbox{where}\ g\in\{B^\prime, w, B^\prime w,\ldots\}. \end{equation}

We include in the statistics only legs that lie entirely within the central region $0.2\leq z\leq 0.8$. In this portion of the domain, the horizontally averaged temperature and salinity profiles are linear in all of the simulations that we performed (table 1). Thus, we have excluded from the leg-based statistics the contribution of the boundary layers.

Leg lengths are compared with the characteristic horizontal size of the doubly diffusive blobs moving in the fluid. A useful quantification of such size is given by the following integral scale (Paparella & von Hardenberg Reference Paparella and von Hardenberg2012):

(2.10)\begin{align} l_{\omega} := {\rm \pi} \left\langle \frac{\displaystyle\sum_{n=1}^{\infty}k_n^{{-}1} |\hat{\omega}|^2(z, k_n)}{\displaystyle\sum_{n=1}^{\infty}|\hat{\omega}|^2(z, k_n)}\right\rangle_{{z}}, \end{align}

where $k_n=2{\rm \pi} n/L_x$ is the horizontal wavenumber, $|\hat {\omega }|^2$ is the power spectral density of the vorticity field (where it is intended that the Fourier transform is only taken along the horizontal direction), and the angular bracket denotes a time and vertical average in the interval $0.2\leq z\leq 0.8$.

3.1. Link between leg length and buoyancy anomaly, vertical velocity, buoyancy flux

Figure 1 shows the bivariate frequency of leg length $\zeta$ and of leg buoyancy anomaly $B^\prime _l$, normalized with respect to the number of legs, for the simulations with $R_S=10^{12}$ and $R_S=10^{14}$. It also shows the conditional expectation $E(B^\prime _l\,|\,\zeta )$, namely the average buoyancy anomaly of legs having length $\zeta$. (Similar figures for the remaining simulations are presented in the supplementary materials.)

Figure 1. Bivariate frequency (normalized with respect to the number of legs, see table 2) of leg buoyancy anomaly ($B^\prime _l$) and leg length ($\zeta$) for (a) $R_S=10^{12}$, (b) $R_S= 10^{14}$. The red curve shows the conditional expectation $E(B^\prime _l\,|\,\zeta )$.

It is immediately evident that long legs do not occur unless their buoyancy anomaly is close to an optimum value (approximately ${\pm }0.8$ for $R_S=10^{12}$ and ${\pm }0.35$ for $R_S=10^{14}$). Short legs, instead, may carry a wide range of buoyancy anomalies. Extreme buoyancy anomaly values are always associated with short legs, but these extreme events are rare: the conditional expectation $E(B^\prime _l\,|\,\zeta )$ reveals that short legs are, on average, associated with small buoyancy anomalies. The sign of the buoyancy anomaly generally agrees with that of the leg length. Thus, as expected, Lagrangian particles associated with positive buoyancy anomaly move upward, and those associated with negative buoyancy anomaly move downward. Given the overall stratification of a fingering convection set-up (heavier fluid underneath lighter fluid), this is the signature of the up-gradient density transport, which is the defining characteristic of doubly diffusive convection. However, short legs make an exception: for legs having length closer to zero, the sign of the buoyancy anomaly is, on average, the opposite of that of the leg length (see the wiggle in $E(B^\prime _l\,|\,\zeta )$ near $\zeta =0$), suggesting that short legs, on average, are associated with a down-gradient flux.

Figure 2 shows the bivariate frequency of leg length $\zeta$ and of leg vertical velocity $w_l$. Just as for buoyancy anomaly, long legs are associated with a relatively narrow band of vertical velocities (either positive, for up-going legs, or negative, for down-going ones). Short legs, instead, exhibit a wider range of vertical velocities, including velocities exceeding those of the long legs. However, on average, a shorter leg will move slower.

Figure 2. Bivariate frequency (normalized with respect to the number of legs, see table 2) of leg vertical velocity ($w_l$) and leg length ($\zeta$) for (a) $R_S=10^{12}$, (b) $R_S= 10^{14}$. The red curve shows the conditional expectation $E(w_l\,|\,\zeta )$.

Figure 3 shows the bivariate frequency of leg length $\zeta$ and of leg advective buoyancy flux $(B^\prime w)_l$. Once again, we find that long legs are associated with a relatively narrow band of advective buoyancy fluxes. Furthermore, long legs always carry positive (up-gradient) fluxes, whereas short legs may carry both positive and negative (down-gradient) fluxes.

Figure 3. Bivariate frequency (normalized with respect to the number of legs, see table 2) of leg advective buoyancy flux ($(B^\prime w)_l$) and leg length ($\zeta$) for (a) $R_S=10^{12}$, (b) $R_S= 10^{14}$. The red curve shows the conditional expectation $E((B^\prime w)_l\,|\,\zeta )$.

Figures 4 and 5 show all three bivariate distributions ($\zeta$ versus $B^\prime _l$, $w_l$ and $(B^\prime w)_l$) for respectively the simulations with $R_\rho =2.0$ and $R_\rho =2.8$, and $R_S=10^{12}$ (those for the other density ratios are in the supplementary materials). Even at high density ratios, the Lagrangian dynamics appears to show the same qualitative properties as in the runs at $R_\rho =1.2$. At $R_\rho =2.0$, and for long legs, the fluctuations around the conditional expectation $E(B^\prime _l\,|\,\zeta )$ are the smallest compared with any other parameter set that we examined. However, at $R_\rho =2.8$, the spread around the average is much larger.

Figure 4. Bivariate frequency (normalized with respect to the number of legs, see table 3) of (a) leg buoyancy anomaly ($B^\prime _l$), (b) leg vertical velocity ($w_l$), (c) leg advective buoyancy flux ($(B^\prime w)_l$) and leg length ($\zeta$) for $R_S=10^{12}$, $R_{\rho }=2.0$ . The red curve shows the conditional expectations $E(B^\prime _l\,|\,\zeta )$, $E(w_l\,|\,\zeta )$, $E((B^\prime w)_l\,|\,\zeta )$.

Figure 5. As in figure 4, but for $R_S=10^{12}$, $R_{\rho }=2.8$.

We approximate the conditional expected flux at a given leg with the following expression:

(3.1)\begin{equation} E\left((B^\prime w)_l\,|\,|\zeta|\right) \approx F_{max}\left(1-\exp\left({-\frac{|\zeta|}{\zeta_0}}\right)\right), \end{equation}

where the free parameters $\zeta _0$ and $F_{max}$ are least-squares fitted to the simulation data. The value of $F_{max}$ estimates the flux that would be associated with legs of formally infinite length. In figure 6, we compare the conditional flux $E((B^\prime w)_l\,|\,|\zeta |)$ of all the simulations of table 1 by normalizing it with respect to $F_{max}$, and plotting it as a function of $|\zeta |/l_\omega$, which is the unsigned leg length expressed in units of the blob size of each simulation. The results of the simulations with $R_\rho =1.2$ collapse almost perfectly on top of each other. The flux reaches half of its maximum for leg lengths of approximately 10 times the size of a blob, and saturates at $|\zeta |/l_\omega \approx 40$. For $R_S=10^{12}$, the curves become progressively less steep as the density ratio $R_\rho$ is increased, and at high density ratio, saturation is not attained. The simulation with the highest density ratio $R_\rho =2.8$, however, breaks the trend, and saturates more quickly than any other.

Figure 6. Normalized conditional expectation $E((B^\prime w)_l\,|\,|\zeta |)/F_{max}$ as a function of the leg length in units of the blob size $|\zeta |/l_\omega$ for all the simulations in table 1. Curves labelled with the density ratio $R_\rho$ refer to simulations with $R_S=10^{12}$. Curves labelled with the haline Rayleigh number $R_S$ refer to simulations with $R_\rho =1.2$.

3.2. What carries the flux

The bivariate distributions presented above show that all long legs tend to be associated with roughly the same values of buoyancy anomaly ($B_l^\prime$), vertical velocity ($w_l$) and up-gradient advective buoyancy flux ($(B^\prime w)_l$). However, they do not immediately reveal how the flux is partitioned among legs of different lengths. In other words, we still need to ascertain whether most of the flux in fingering convection is carried by portions of the fluid that coherently travel in the same direction, or, instead, the flux is mostly determined by the jittery, incoherent motion revealed by the large amount of short legs. To answer this question, we introduce the cumulative flux fraction, $C_F(|\zeta |)$, which is the fraction of the total flux associated with legs having length up to $|\zeta |$:

(3.2)\begin{equation} C_F\left(|\zeta|\right) := \frac{\displaystyle\int_0^{|\zeta|} \int_{-\infty}^{+\infty} p(\xi, F) F\,\mathrm{d}F\,\mathrm{d}\xi}{\displaystyle \int_{0}^{\infty}\int_{-\infty}^{+\infty} p(|\zeta|, F) F\,\mathrm{d}F\,\mathrm{d}|\zeta|}, \end{equation}

where $p$ is the joint probability distribution of unsigned leg length $|\zeta |$ and leg advective buoyancy flux $F=(B^\prime w)_l$.

The cumulative flux fraction is shown in figure 7 for the simulations with $R_S=10^{12}$ and $R_S=10^{14}$ (see supplementary materials for the others). At $R_S=10^{12}$, the cumulative flux is negative (down-gradient) up to lengths $\sim$1.6 times the blob size $l_\omega$. At $R_S=10^{14}$, this value increases to $\sim$2.5 $l_\omega$. The very abundant short legs carry a down-gradient flux which reaches 12 % of the total flux at $R_S=10^{12}$ and 21 % of it at $R_S=10^{14}$ (minima of the cumulative flux fraction). At $R_S=10^{12}$, 50 % of the up-gradient flux is carried by legs longer than $\sim$4.5 $l_\omega$, and this increases to 60 % at $R_S=10^{14}$. Thus, the bulk of the up-gradient flux is carried by legs substantially longer than the size of a blob. Even though these long legs are fewer in number (the distribution of leg lengths is shown in red in figure 7), their positive flux is sufficient to overcome the negative one carried by the much more abundant short legs.

Figure 7. Cumulative flux fraction (green line) as a function of the leg length (left axis). Probability density function of the leg length (red line, right axis, log scale). (a) Simulation with $R_S=10^{12}$, (b) simulation with $R_S=10^{14}$.

In figure 8, we show the cumulative flux fraction and the leg length probability density function for the simulations with $R_S=10^{12}$ and $R_\rho =1.2,\,2.0,\,2.8$. In this case, on the horizontal axis, we show the leg length as a fraction of the domain height. As the density ratio increases, the down-gradient flux carried by short legs becomes negligibly small. This result suggests that down-gradient fluxes are the result of mechanical mixing operated by fast-moving blobs: vertical velocities dramatically decline (figures 2(a), 4(b), 5(b)) with increasing density ratio, and with them, the ability of the blobs to stir the surrounding fluid. Correspondingly, long legs become longer and more abundant (table 3) because blobs are less disturbed by nearby ones as they move vertically. However, any trend associated with lengthening leg lengths breaks down when the density ratio becomes close enough to the marginal stability threshold $R_\rho =Le$ (recall that here we are using $Le=3$). At $R_\rho =2.8$, the exponential tail of the leg length distribution becomes steeper than at $R_\rho =1.2$; very long legs become rarer (table 3); the normalized conditional expectation $E((B^\prime w)_l\,|\,|\zeta |) / F_{max}$ becomes steeper (figure 6). Considering that at marginal stability, the fluid would be at rest, and legs would have zero length, the breakdown of the trends appears to be a natural consequence of the continuous dependence of the flow statistics on the density ratio.

Figure 8. As in figure 7(a), but for the simulations with $R_S=10^{12}$ and density ratio $R_\rho =1.2$ (solid lines), $R_\rho =2.0$ (dashed lines), and $R_\rho =2.8$ (dotted lines). Note that the horizontal axis now expresses the leg lengths as a fraction of the domain height, rather than in blob size units.