4.1. Flow structures of interleaving

We first show the detailed flow structures for case 2 in figures 5(a) and 5(b) by the three-dimensional volume rendering of the streamwise and vertical velocities. A stack of interleaving layers is distinctively observed in the domain. For this case seven layers are visible. The layers extend upward towards the positive streamwise direction and have alternating flow directions with a positive or negative streamwise velocity, which is a clear indication of interleaving. The layers have long horizontal length and relatively small vertical thickness. In figure 5(c) we plot the profiles of the temperature and salinity anomalies at seven different streamwise locations. Here, the anomaly of scalar $\zeta$ is calculated as $\zeta '=\bar {\zeta }-\overline {\bar {\zeta }}$. Hereafter, one bar stands for the spanwise average, and two bars for the spanwise and temporal average. Coherent temperature and salinity anomalies exist along the interleaving currents.

Figure 5.

Thermohaline interleaving for case 2 at $t=35\,000$ depicted by the volume renderings of (a) the streamwise velocity and (b) the vertical velocity, in which the positive $y$-direction is from left to right. (c) The profiles from left to right are for temperature anomaly (red dashed line) and salinity anomaly (blue solid lines) at $y=iL/8$, with $i=1\cdots 7$.

Once the flow reaches the steady quasi-periodic state, the mean temperature and salinity fields drive different types of flows within the domain. To demonstrate this, we plot in figure 6(a) the contours of spanwise and temporally averaged temperature and salinity fields. At the upper left and upper right corners, the isotherms and isohalines are almost horizontal and parallel. The gradients of both scalars are mainly in the vertical direction, and DDC occurs with the diffusive type near the left corner and the fingering type near the right corner. The fingering DDC generates the small-scale vertical motions on the right in figure 5(b). At the left part of the bulk region, where the layering is the strongest, the isotherms and isohalines have different inclination angles. By averaging over the main intrusion region, as marked by the dashed rectangular in figure 6(a), the mean inclination angle for isotherms is approximately $25^{\circ }$ and that for isohalines is approximately $11^{\circ }$. The interleaving currents have an inclination angle of approximately $3^{\circ }$, which is smaller than those of the isotherms and isohalines, a phenomenon consistent with the intrusion theory (Ruddick Reference Ruddick2003; Radko Reference Radko2013).

Figure 6.

(a) Contours of spanwise and temporally averaged temperature $\overline {\bar {T}}$ (solid lines) and salinity $\overline {\bar {S}}$ (dashed lines). (b) Flow regions susceptible to different types of double diffusion. Colours represent the Turner angle $Tu$, with blue colours indicating the diffusive type and red colours the fingering type, respectively. Here, $Tu$ is calculated by the spanwise averaged scalar gradients $\partial _z\bar {T}$ and $\partial _z\bar {S}$ for the flow field shown in figure 5.

The different types of DDC can be also identified by the Turner angle, which is defined by the vertical gradients of the mean scalars as $Tu=135^{\circ }-\arg (\beta _S \partial _z \bar {S} + i\beta _T \partial _z \bar {T})$ (Ruddick Reference Ruddick1983). Here, $i$ is the imaginary unit, and $\arg$ stands for the argument of the imaginary number. Fingering DDC usually happens in the range $45^{\circ }< Tu<90^{\circ }$, and diffusive DDC in the range $270^{\circ }< Tu<315^{\circ }$. In figure 6(b) we present the distribution of $Tu$ for the flow field shown in figure 5, where the mean gradients are calculated after spanwise averaging. Again, both fingering and diffusive DDC exist in the flow field. The region with $45^{\circ }< Tu<90^{\circ }$ agrees with the vertical motions in figure 5(b), namely, those vertical motions are mainly fingering convection. At the upper-right corner, the fingering DDC is driven by the vertical salinity and temperature gradients imposed by the top boundary. In the left part of the bulk, the layers with fingering DDC correspond to the interfaces below the leftward moving currents with higher temperature and salinity and above the rightward moving ones with lower temperature and salinity. Meanwhile, diffusive DDC dominates the upper left region, where $270^{\circ }< Tu<315^{\circ }$. The layers with $270^{\circ }< Tu<315^{\circ }$ correspond to the interfaces above the leftward moving warm and salty currents and below the rightward moving cold and fresh ones. Different types of DDC between interleaving currents are characteristic phenomena of thermohaline intrusion (Ruddick Reference Ruddick2003; Radko Reference Radko2013).

4.2. The upward shift of the layers

Moreover, our numerical results reveal that the vertical location of each layer is not constant. Rather, these layers shift upward over a very large time scale, which is clearly shown in supplementary movie 1 available at https://doi.org/10.1017/jfm.2021.753. Figure 7 displays the time evolution of the vertical profiles of streamwise velocity at location $y=3/8L$, i.e. where the interleaving currents are strongest. For case 2, it takes approximately $10\,000$ non-dimensional time units for one layer to move from the bottom to the top surface. Detailed investigation indicates that the total density difference between the two ends in the streamwise direction periodically changes sign, and close to the bottom boundary new layers are generated accordingly with alternating movement directions.

Figure 7.

The upward shifting of the interleaving layers depicted by the time evolution of the spanwise-averaged profiles of the streamwise velocity at the streamwise location $y=3L/8$ for case 2.

From the contours in figure 7 we extract the upward-shifting velocity $U_s$ of the layers. Specifically, within the height range $0.2< z/H<0.8$, we calculate the auto-correlation coefficient of $\bar {v}$ for the height separation $\delta z$ and the time separation $\delta t$. Then, the shift velocity $U_s$ is determined by the location of the maximal auto-correlation coefficient on the $\delta z$–$\delta t$ plane. The values are given in table 2. Compared with the current velocity, the shift velocity $U_s$ is much smaller, by approximately two orders of magnitude.

The upward shift is caused by the continuously emerging new layers near the bottom boundary, which are generated by the periodically oscillating density difference between the two ends of the domain. This lateral asymmetry may be the key reason for the upward migration. It should be mentioned that the linear instability analysis revealed a migration of the unstable modes in the heated sidewall configuration (Kerr Reference Kerr1989). We calculate the difference of mean density over two streamwise locations $y_1=7L/8$ and $y_2=L/8$, namely $\delta \rho =\langle \rho '\rangle _{y_1}-\langle \rho '\rangle _{y_2}$. Here, $\langle \,\rangle _y$ stands for the spatial average over the streamwise cross-section at the location $y$. The density is calculated as $\rho '=\varLambda S - T$. The total density flux through the top surface is measured accordingly as $F_\rho ^{top}=\kappa _T \langle \partial _z T\rangle _{top} - \kappa _S \varLambda \langle \partial _z S\rangle _{top}$. In figure 8(a) we plot the time histories of $\delta \rho$ and $F_\rho ^{top}$ for case 2. Both quantities fluctuate with a similar period. Moreover, the variation of $\delta \rho$ lags behind $F_\rho ^{top}$ by a phase difference of approximately ${\rm \pi} /2$. From the time history of $\delta \rho$ we can extract the mean period of oscillation $t_{density}$. If we assume that for each period two layers develop near the bottom boundary, then the upward-shifting velocity $U_s$ should be very close to the velocity scale $2h/t_{density}$. This is indeed the case, as shown in figure 8, where all data points are very close to the line $U_s=2h/t_{density}$.

Figure 8.

(a) The time histories of density flux through the top surface $Nu_\rho ^{top}$ (black dash-dotted line and left ordinate) and the mean density difference $\delta \rho$ (red solid line and right ordinate) between $y=7L/8$ and $L/8$ for case 2. (b) The upward-shifting velocity $U_s$ vs $2h/t_{density}$ with $h$ being the layer thickness and $t_{density}$ the period of the $\delta \rho$.

4.3. Statistics of interleaving layers

One of the key parameters of intrusion is the layer thickness $h$. Previous experiments by Ruddick & Turner (Reference Ruddick and Turner1979) (RT experiments) found the scaling law of $h_{RT}$, which is proportional to the so-called Chen scale (Chen Reference Chen1974) developed for the bounded intrusion model. Later, Toole & Georgi (Reference Toole and Georgi1981) obtained another scale $h_{TG}$ for the unbounded case. These two scales were then proven to be the two end limits in asymptotic analysis for the intrusion front width, which was demonstrated by Niino (Reference Niino1986) and Simeonov & Stern (Reference Simeonov and Stern2004) in different ways. The latter also introduced an extra scale $h_{SS}$ for an intermediate state. Here, we will first use Reference NiinoNiino's (Reference Niino1986) method to reveal that the fronts are narrow enough in our cases (the same as in the RT experiments) and then demonstrate that the layer thickness in our simulations indeed has a scaling law which is essentially the same as that for $h_{RT}$. It should be pointed out that the model developed in Simeonov & Stern (Reference Simeonov and Stern2004) requires a preset vertical density ratio $\varLambda _z$ and its relation with Nusselt number $Nu(\varLambda _z)$ for the fingering regime. In our configuration, however, both fingering and diffusive regimes are present simultaneously, so it is difficult to estimate $\varLambda _z$ or $Nu(\varLambda _z)$. Actually, the Niino (Reference Niino1986) theory also applies to the fingering regime, but he demonstrated that the theory is valid for the diffusive cases as far as the salt fingers are present, which is satisfied in our simulations.

The following parameters are required in the Niino (Reference Niino1986) theory: the frontal stability parameter $G=(g(1-\gamma _f) \beta _S {\rm \Delta} ^{h}_S/2)^{6}/(K_S^{2} l^{2} N^{10})$, the modified Schmidt number $\epsilon =\nu /K_S$, the stratification parameter $\mu =g(1-\gamma _f)\beta _S \partial _z\overline {\bar {S}}/N^{2}$ and the modified stability parameter $R=G/\epsilon$. Here, $\gamma _f$ is the flux ratio for fingering DDC and we use the constant value $0.88$ since the diffusive ratio $\tau$ equals $1/3$, which was also used in RT experiments. Also, $K_S$ is the eddy diffusivity of salinity and we calculated it in the same way as $Nu_S$ (see § 4.4); $l$ is the half-width of the intrusion front and we assume it to be of the same order as $L/2$ since the interleaving layers nearly occupy the whole bulk area in all cases; ${\rm \Delta} ^{h}_S$ is the lateral variation of salinity within the water mass; and $N=\sqrt {-(g/\rho )\partial _z\rho }$ is the buoyancy frequency. In the RT experiments, both the lateral salinity variation and the vertical scalar stratification are carefully set as the control parameters, which is not the case in our simulations. Instead, we take ${\rm \Delta} ^{h}_S$ as the averaged salinity difference between two endwalls in the streamwise direction, and calculate $N$ and $\partial _z\overline {\bar {S}}$ by the mean density and salinity difference between the top and bottom boundaries. The results are summarized in table 3. The fronts can be considered narrow when $R<40(1+\mu )^{5.4}$, and obviously our results fall into this range. Thus, $h$ should have the same scaling law for $h_{RT}$. Note that $\mu <0$ indicates the diffusive regime in our cases. As mentioned above, Niino (Reference Niino1986) demonstrated that the theory is still valid for $-1<\mu <0$ as long as salt fingers are present. The values of $G$ and $R$ are much smaller than those in the RT experiments (Niino Reference Niino1986). This may be attributed to the small ${\rm \Delta} ^{h}_S$, which is qualitatively consistent with the Simeonov & Stern (Reference Simeonov and Stern2004) theory.

Table 3.

Summary of the key parameters in the Niino (Reference Niino1986) model. Columns from left to right are as follows: case number, the lateral variation of salinity across the front, the frontal stability parameter, the modified Schmidt number, the stratification parameter, the modified stability parameter and the low and high critical values for determination of the front width. When $R<40(1+\mu )^{5.4}$, the front can be considered narrow and the layer thickness has a scaling law the same as $h_{RT}$ (Ruddick & Turner Reference Ruddick and Turner1979). When $R>2\times 10^{5}(1+\mu )^{4.9}$, the front can be considered wide and the layer thickness has a scaling law the same as $h_{TG}$ (Toole & Georgi Reference Toole and Georgi1981).

Now we would like to examine the scaling law of layer thickness, which is sampled from individual spanwise mean profiles at $y=3L/8$ and averaged over time and all layers. Since in the present flows the fronts are narrow and the layer thickness should follow the RT scaling $h_{RT} = C_h(g\beta _S{\rm \Delta} ^{h}_S)/N^{2}$ (Ruddick & Turner Reference Ruddick and Turner1979; Niino Reference Niino1986), in figure 9(a) we plot $h$ vs $(g\beta _S{\rm \Delta} ^{h}_S)/N^{2}$. Indeed the data follow the scaling law with $C_h\approx 0.066$.

Figure 9.

Properties of the interleaving layers. (a) Layer thickness $h$ vs $(g\beta _S{\rm \Delta} ^{h}_S)/N^{2}$. (b) The maximal current velocity $v_m$ vs $Nh$. (c) The horizontal heat flux $Nu_T^{h}$ vs $a^{-1/2}$ with $a$ being the isohaline slope.

Another parameter that has been widely studied is the maximum horizontal velocity $v_m$ of the interleaving layers. For the different frontal widths and DDC types, previous experiments (Ruddick & Turner Reference Ruddick and Turner1979) and simulations (Simeonov & Stern Reference Simeonov and Stern2007,  Reference Simeonov and Stern2008) all observed the same scaling law $v_m\sim C_v N h$ with different values for the prefactor $C_v$. We use the spanwise-averaged profiles at $y=3L/8$ to find the maximal horizontal velocity $v_m$ of all layers during the whole sampling period. The result is consistent with the previous work and the prefactor $C_v$ is approximately $0.15$, as shown in figure 9(b). The value of $C_v$ is very close to those reported in Simeonov & Stern (Reference Simeonov and Stern2007,  Reference Simeonov and Stern2008).

As shown in figure 6(a), the layering intrusion in the present flow is of diffusive type. For such a type of intrusion, the horizontal Nusselt number of temperature $Nu_T^{h}$ can be related to the isohaline slope $a$ by the scaling law $Nu_T^{h} \sim a^{-1/2}$ (Simeonov & Stern Reference Simeonov and Stern2008). Here, $Nu_T^{h}$ represents the ratio of the horizontal heat flux to the horizontal temperature gradient, namely $Nu_T^{h}=\langle vT \rangle _h/( \kappa _T{\rm \Delta} _T^{h}{L}^{-1} )$. The average $\langle \rangle _h$ is conducted over the streamwise cross-section at $y=3L/8$ and over time, and ${\rm \Delta} _T^{h}$ is again taken as the averaged temperature difference between two endwalls in the streamwise direction. The isohaline slope is estimated by averaging over the main intrusion region, as explained in § 4.1. Figure 9(c) displays the dependence of $Nu_T^{h}$ vs $a^{-1/2}$. The data roughly follow the scaling law $Nu_T^{h} \sim 128.8a^{-1/2}$, where the prefactor is fixed by linear fitting.

4.4. The global transport

Finally, we proceed to analyse the global transport properties in the horizontal direction. Due to our specific flow configuration, both heat and salinity are transported into the domain from the right part of the top surface with higher temperature and salinity, then move horizontally towards the left side within the domain bulk and finally exit the domain over the left part of the top surface with lower temperature and salinity. Such behaviours are shown in figure 10 by the distributions of local Nusselt number over the top boundary for case 2. The local Nusselt number is calculated as $Nu^{l}_\zeta = (-\kappa _\zeta \partial _z \overline {\bar {\zeta }})/( \kappa _\zeta {\rm \Delta} _\zeta {L}^{-1})$ with $\zeta =T$ or $S$. The zero point of $Nu^{l}$ is not at the middle point $z=L/2$ but shifted towards the warm salty end, which happens for all the four cases considered. Near the left end the vertical flux is generated by the diffusive DDC process, while near the right end it is generated by the fingering DDC. In between, the positive and negative peaks next to the zero point are caused by the upward-shifting layers when they reach the top surface.

Figure 10.

The local Nusselt number $Nu^{l}$ over the top boundary for case 2.

Moreover, the total flux through the entire top surface oscillates around zero with a very long period. For example, case 1 has a time period of approximately 3000 time units, as shown in figure 2. Nevertheless, the temporally averaged total flux should be approximately zero when the system is in the statistically steady state and the average is calculated over a long enough time period. We measure these global fluxes $Nu_\zeta$ by integrating $Nu^{l}_\zeta$ over the range at either side of the zero-flux point. In addition, we define the horizontal Rayleigh number $Ra_L=(g\beta _T L^{3} {\rm \Delta} _T)/(\nu \kappa _T)$ to account for the effect of the streamwise length $L$. The dependences of $Nu_T$ and $Nu_S$ on $Ra_L$ are plotted in figures 11(a) and 11(b), respectively. Both Nusselt numbers follow a very similar scaling law: $Nu_\zeta \sim Ra_L^{\alpha _\zeta }$ with $\alpha _\zeta$ very close to $0.2$. The Reynolds number $Re=U_{rms}L/\nu$ follows a scaling law of $Re \sim Ra_L^{0.44}$. Here, $U_{rms}$ is the root-mean-square value of the velocity magnitude.

Figure 11.

The global responses vs the horizontal Rayleigh number $Ra_L$. (a) The global heat flux. (b) The global salinity flux. (c) The Reynolds number. (d) The horizontal density flux ratio $\gamma$. In (a,b) the error bars indicate the maximal deviation from zero of the total fluxes through the top surface.

Since both the heat and salinity fluxes have very similar scaling behaviours, the density flux ratio $\gamma =(\beta _T\langle v T\rangle _0)/(\beta _S\langle v S\rangle _0)$ should be nearly constant for different $Ra_L$. The average $\langle \,\rangle _0$ is conducted over the streamwise cross-section at the zero-flux point of the top boundary, and over time. Figure 11(d) shows that $\gamma$ lies between $1.7$ and $1.8$. Nearly constant $\gamma$ implies that the ratio of the density anomaly flux due to the heat flux to that due to the salinity flux exhibits very weak dependence on $Ra_L$. Moreover, $\gamma >1$ indicates that the net density anomaly flux is in the same direction as the scalar gradients along the top boundary.