3.1. Definitions and metrics

We will employ the notation of § 2.2, as sketched in figure 2, with dimensionless quantities based on the flume scales (2.4a–c). Here, the bottom wall of the flume is located at $y=0$, and $y=\eta$ denotes the top of the bed. However, at times it will be more meaningful to compare dimensionless quantities based on the free depth scales, which can be obtained by normalizing the dimensionless quantities based on the flume scales with $T_H$ and $U_H$, as defined in (2.6) and (2.5), respectively. The origin is defined in space as the position of the lock gate in and in time as the original lock lift time. The experiments are only analysed while they are in the slumping phase, i.e. before the lock runs out of fluid to constantly drive the currents.

We define the buoyant height $h$ by the integral over the concentration above the bed (Shin et al. Reference Shin, Dalziel and Linden2004)

(3.1)\begin{equation} h(x,t) = \int_{\eta}^1 \langle c \rangle_z\, \textrm{d}y, \end{equation}

where the angled brackets denote the average in the direction of the subscript, i.e. $\langle c \rangle _z$ is the concentration averaged in the spanwise $z$-direction. This height is indicative of the buoyant driving force within the current, depending on both the salinity within the current and the vertical extent of the current. The front location $x_F(t)$ of the gravity current is then determined as the rightmost location at which $h(x,t) = 0.1$, so that we subsequently obtain the front speed $v_F(t)$ as

(3.2)\begin{equation} v_F(t)= \frac{\textrm{d} x_F(t)}{\textrm{d}t} . \end{equation}

We furthermore define the current envelope $h_c(x,t)$ via the maximum $y$-locations where the spanwise-averaged density has the value 0.02. We observed repeat experiments utilizing single camera windows to identify the errors associated with these parameters. It was found for trends in these parameters the range of these repeats did not exceed 3 %.

3.2. Overview of the gravity current development

The $c=0.1$ salinity contours shown in figure 3 demonstrate the spatio-temporal development of the flow for the representative simulation S3L15B, with $d_p=0.066$ and $Re_H = 10\ 355$. At $t/T_H =3.3$, the current displays coherent spanwise Kelvin–Helmholtz vortices, along with regular, small frontal lobes that move along the grooves between the top layer of spheres and hence reflect the periodicity of the bed arrangement. At $t/T_H=7.7$ these lobes remain nearly periodic, although their wavelength has doubled, while the tail is already fully three-dimensional. By $t/T_H=16.5$ the entire current has transitioned to a turbulent state.

Figure 3.

The $c=0.1$ salinity isosurface for simulation S3L15B at times $t/T_H=3.3$, 7.7 and 16.5. While the early stage of the current is dominated by coherent spanwise vortices, the flow transitions to a fully turbulent state later on.

Figure 4(a) displays the instantaneous salinity field in the plane $z=0.4825$ for $t/T_H=17.6$, along with the projected velocity vectors. The current is seen to mix vigorously with the ambient fluid along the upper interface. The close-ups in (b,c) focus on the detailed flow features near the top of the particle bed, and indicate the numerical resolution of 40 grid cells per particle diameter, by showing the velocity vectors and salinity values at every grid point. The close-ups furthermore visualize the fluid exchange between the current and the bed, as dense current fluid leaks through the gaps between the spheres into the pore spaces. In addition, they show that the flow separates from the top of the spheres, which results in the strong deformation and folding of the fluid interface and thereby promotes mixing.

Figure 4.

Concentration $c$ (in blue) and velocity vectors in the plane $z=0.4825$, for simulation S3L15B at $t/T_H=17.6$. Panel (a) demonstrates that the current mixes vigorously with the ambient fluid along the upper interface, while it also exchanges fluid with the porous bed. The close-ups in (b,c), indicated by black rectangles in (a), highlight details of this current–bed interaction. Panel (b) shows every grid cell, so that it provides a measure of the spatial resolution.

Figure 5 displays the salinity concentration field averaged in the spanwise $z$-direction over the fluid region only, $\langle c \rangle _{z}/ \langle \xi _f \rangle _z$, for simulation S3L15B at times $t/T_H=$ 7.7, 12.1 and 26.4. The current front propagating above the particle bed is qualitatively similar to that advancing over a flat bottom (Härtel et al. Reference Härtel, Meiburg and Necker2000). In addition, we observe a second front inside the porous bed that moves much more slowly, which is consistent with earlier experimental observations by Cenedese et al. (Reference Cenedese, Nokes and Hyatt2018). This second, slower front is driven by the dense fluid that leaks into the bed along the entire length of the current. The small difference between the salinity concentration values above and inside the bed in the lock region is due to numerical errors associated with the volume-of-fluid method, which cause a slight diffusive salinity flux into the particles.

Figure 5.

Spanwise-averaged salinity $\langle c \rangle _{z}/ \langle \xi _f \rangle _z$ in the head region of simulation S3L15B, at times $t/T_H$= 7.7 (a), 12.1 (b) and 26.4 (c). A conventional gravity current front propagates above the bed, while a second front advances more slowly inside the bed. With time, both the buoyant height $h$ (long–short dashed line) and the envelope $h_c$ (dashed line) of the current become more uniform along the length of the current.

During the early stages of the flow, the spanwise coherent Kelvin–Helmholtz vortices mentioned earlier mix the two fluids across nearly the entire free depth of the flow, cf. figure 5(a). Later on, after the flow has become fully three-dimensional and turbulent, the most intense mixing is restricted to a narrower region around the interface. Visual inspection of the $c=0.1$ concentration isosurfaces at various intermediate times (not shown) indicates that the simulation with $Re_H=10\ 355$ (case S3L15C) transitions to turbulence somewhat earlier than the one for $Re_H=7486.2$ (case S1L15A), in line with our expectation that higher Reynolds numbers should promote more vigorous turbulence. We note that the dense current fluid quickly enters the spaces between the upper half of the top layer of spheres, as the resident light fluid in these exposed pores is easily carried away by the current. The infiltration of the current fluid into the lower layers of the sediment bed takes much longer, however, as the dense current fluid has to displace the lighter, resident fluid via a Rayleigh–Taylor-like instability in an effective porous medium. With time, both the buoyant height $h$ (long–short dashed line) and the current envelope $h_c$ (dashed line) tend to become more uniform along the length of the current, see figure 5(c) for $t/T_H=26.4$.

Figure 6 displays the salinity concentration field averaged in the spanwise $z$-direction over the fluid region only, $\langle c \rangle _{z}/ \langle \xi _f \rangle _z$, for the experiment E3L15 at times $t/T_H= 7.7$, 12.1, and 26.4. The current front propagating above the particle bed is qualitatively similar to the simulations shown in figure 5. The simulations show a less uniform buoyant height and current envelope at earlier times because the experiments transitioned to a fully turbulent state at earlier times.

Figure 6.

Spanwise-averaged salinity $\langle c^+ \rangle _{z}$ in the head region for experiment E3L15, at times $t/T_H$= 7.7 (a), 12.1 (b) and 26.4 (c). Current propagates above bed similarly to the simulation counterparts in figure 5; however, with a more uniform buoyant height $h$ (long–short dashed line) and the envelope $h_c$ (dashed line).

Occasionally, it is preferable to analyse the flow in the reference frame moving with the front, with the front location at the origin of the coordinate system

(3.3)\begin{gather} c^+( {\boldsymbol{x} } - \boldsymbol{e}_x x_F(t ),t) = c ( {\boldsymbol{x}}, t ), \end{gather}

(3.4)\begin{gather}{\boldsymbol{u}}^+ ( {\boldsymbol{x}} - \boldsymbol{e}_x x_F(t) ,t ) = \boldsymbol{u} ( {\boldsymbol{x} } , t ) - v_F(t) \boldsymbol{e}_x . \end{gather}

We denote variables in this reference frame by a $^+$-superscript. Figure 7 shows the spanwise-averaged concentration $\langle c^+ \rangle _{z}$ for experiment 3L15 at two different times. In addition to averaging in the spanwise direction, which is inherent to the LA method, we average the distribution around the given times by $\pm 1 T_H$ in order to smooth out fluctuations. For comparison, we display the corresponding simulation results for $t/T_H=22$ in (c).

Figure 7.

Spanwise-averaged salinity $\langle c \rangle _{z}$ in the head region for experiment E3L15, at times $t/T_H=5.3$ (a) and 22 (b). The panels also include the buoyant height $h$ (long–short dashed line) and the current envelope $h_c$ (dashed line). Corresponding simulation results for $t/T_H= 22$, shown in (c), are seen to agree closely with the experimental values. All distributions are averaged around the given time for $\pm 1 T_H$.

We note that the early experimental salinity field at $t/T_H=5.3$ does not indicate the presence of the strong, spanwise coherent Kelvin–Helmholtz vortices that we had observed during the initial simulation stages of the simulation, cf. figure 5(a). This is likely a result of the much stronger initial perturbation introduced in the experiment by the removal of the gate, which causes the flow to transition to a fully three-dimensional turbulent state more rapidly than in the simulation. Nevertheless, by $t/T_H=22$ the experimental and computational concentration fields agree closely, especially with regard to the buoyant height as a function of the streamwise coordinate.

3.3. Front position and velocity

A key question concerns the influence of the porous bed on the front velocity of the current. Figure 8 presents experimental results for the non-dimensional front location $(x_F -x_0)/H$ as function of time, where the value of zero corresponds to the front being at the position of the lock. While (a) shows values for the single-layer flows and the flat bottom, (b) provides corresponding data for the three-layer flows. The two-layer case behaves qualitatively similar, so that it is not shown. Within each panel the water depth $\tilde L_y$ is varied in order to assess its influence. Table 1 lists all experimental parameters.

Figure 8.

Front location and velocity as functions of time, for the experiments listed in table 1. The left (right) column shows results for a single layer (three layers) of spheres. The upper (lower) row presents the front locations (front velocities) in free depth scales.

The experimental time $\tilde t$ was normalized with the free depth time $\mathcal {T}_H$, and the front position with the free depth $\tilde H$. The time $\tilde t=0$ corresponds to the moment when the lock is opened. For a constant number of layers, this scaling causes the curves for the front locations to collapse approximately, although the shallower cases still propagate somewhat more slowly even when scaled with these free depth units, reflecting their effectively lower Reynolds numbers. In particular, we note that the current advancing over a smooth bottom travels only a few percentage points faster than the one propagating over a single layer of spheres. This indicates that during the early stages of the flow, the bottom roughness and the loss of current fluid into the porous bed of a single layer of spheres do not significantly retard the current.

Differentiation of the location values in (a,b) yields the front velocities, shown in (c,d), normalized with $\mathcal {U}_H$, as a function of time. For all water depths the front velocity $\tilde v_F/ \mathcal {U}_H$ decreases with time, although this trend is most pronounced for the shallow water depths, which indicates that over time the bed roughness and porosity have a relatively larger influence on shallower currents.

Very prominent features of the front velocity data are their strong oscillations in time, which to the best of our knowledge had not been reported in earlier laboratory experiments. We hypothesize that these front velocity oscillations are caused by the presence of seiche waves generated by the removal of the lock gate. In order to test this hypothesis, let us consider the first seiche mode in a channel, which has a wavelength $\lambda$ of twice the channel length. For our experiments this wavelength is large compared to the water depth, so that we can consider the seiche a shallow, linear wave with a phase speed $\tilde c_w = \sqrt {g\tilde H}$. The period $\tilde T_w$ of the wave, and the maximum speed $\tilde u_w$ of a fluid particle beneath the surface are given by

(3.5)\begin{gather} \tilde T_w = \frac{\tilde \lambda}{\tilde c_w} = \frac{2 \tilde L_x}{\sqrt{g\tilde H}}, \end{gather}

(3.6)\begin{gather}\tilde u_w = \frac{g \tilde \chi_w}{\tilde c_w} = \frac{g\tilde \chi_w}{\sqrt{g \tilde H}} , \end{gather}

where $\tilde \chi _w$ denotes the maximum free surface deflection of the wave, which is small enough to be hardly visible in the experiment. For experiment E3L15 ($\tilde {H} = 0.126$ m and $L_x = 6.2$ m) (3.5) gives a period of $\tilde {T}=11.2$ s, or $\tilde {T}/\mathcal {T}_H = 7.4$. This value closely matches the oscillation period observed for E3L15 in figure 8. For the other experiments we found a similar level of agreement. For E3L15 the maximum velocity was 0.5 cm s$^{-1}$ greater than the mean. We employ this value as an estimate for the maximum velocity of a fluid particle beneath the surface in the absence of the current. Equation (3.6) then yields a maximum free surface deflection of 0.6 mm, which is small enough not to be noticed in most laboratory set-ups. In the following discussion, we assume that the seiche waves did not impact the mean current behaviour, beyond causing the current velocity to oscillate. We furthermore remark that the simulations do not generate seiche modes, as they assume a rigid upper surface.

Figure 8(c,d) includes linear fits to the front velocity data over the entire time $\tilde T_{exp}$ of the experiments in order to better visualize the trends in the data. We fit (by means of minimizing the squared error $\epsilon ^2$) the slope $\alpha$ to the experimental data, which formally reads as

(3.7)\begin{equation} \tilde v_F / \mathcal{U}_H = \langle \tilde v_F \rangle_t/\mathcal{U}_H+ \alpha (\tilde t - \tilde T_{exp}/2 -\tilde t_0)/\mathcal{T}_H + \epsilon.\end{equation}

In (3.7), the average front velocity is obtained as $\langle \tilde v_F \rangle _t = (\tilde x_F(\tilde T_{exp})- \tilde x_F(\tilde t_0)) / \tilde T_{exp}$, where $\tilde t_0$ is the time when the front is detected initially. Moreover, we define the densimetric Froude number

(3.8)\begin{equation} Fr(t)=\tilde v_F(t) / \mathcal{U}_H. \end{equation}

Column 9 of table 1 provides the average dimensionless front speeds. For most experiments the front velocity decreases with time, although for the deeper cases with only one particle layer, we observe a small increase during the analysed time interval. In general, the currents slow down more rapidly for larger relative bed heights $\tilde \eta / \tilde H$.

We perform a corresponding fit for the simulations, although there we exclude the initial laminar phase of the flow. The time interval $t\in \tau =[\tau _0, \tau _1]$ over which we evaluate the linear fit is given in table 3. In order to conduct meaningful comparisons, we evaluate the linear fit (3.7) at identical times $t/T_H=20$ for all experiments and simulations, to obtain $Fr|_{t/T_H=20}$. When plotting $Fr|_{t/T_H=20}$ against $d_p/H$, see figure 9, with error bars included based on the range of repeat experiments, we find that the Froude number is primarily a function of the relative particle size $d_p/H$, whereas the number of layers generally plays only a secondary role. An obvious exception to this rule is the single-layer case (square symbol) for large $d_p/H$, which differs substantially from the two- and three-layer cases with a similar $d_p/H$ (star and diamond symbol), in line with the observations reported in § 3.4.

Figure 9.

Experimental and numerical results showing the correlation between sphere diameter $d_p/H$ and Froude number evaluated at $t/T_H=20$, $Fr|_{t/T_H=20}$. Different symbol shapes represent different numbers of layers, with red (black) symbols representing experiments (simulations). Error bars for experiments are based on the range of repeat experiments. Simulations for constant $d_p/H$ but changing Reynolds and Schmidt numbers are included. The precise parameter values can be found in table 3 (simulations) and table 1 (experiments). As a general trend, the front velocity decreases with increasing particle size, whereas the number of layers plays a secondary role.

Cenedese et al. (Reference Cenedese, Nokes and Hyatt2018) carried out lock-exchange experiments with similar depths to those in the experiments reported here. They found Froude numbers similar to ours for currents propagating over smooth beds. These are also similar to other experimental studies of gravity currents over smooth beds (Simpson Reference Simpson1999; Shin et al. Reference Shin, Dalziel and Linden2004; Cantero et al. Reference Cantero, Balachandar and Garcia2007). Cenedese et al. (Reference Cenedese, Nokes and Hyatt2018) also carried out experiments with staggered beds of vertically arranged circular cylinders and when these were densely distributed the currents propagated atop the roughness. With cylinder heights of 10 mm, the same roughness height as our experiments with 1 layer of particles, they observed Froude numbers consistently lower than ours. We explain this by their experiments having a larger and more unobstructed volume of fluid within the bed such that the vertical exchange of fluid between the bed and the current head was larger than in our experiments.

Figure 10(a) displays the front velocities of experiment E3L15 and the corresponding simulation S3L15C as functions of time.

Figure 10.

Comparison of experimental and numerical results for the front velocities/ densimetric Froude numbers. In (a) the evolution for experiment E3L15 is compared to different simulations showing a close fit between the experiment and simulations. Panel (b) compares experiments E0L20, E3L20, E1L15, E3L15 and E3L10 with simulations S0L15B, S3L20A, S3L15B and S3L10A in terms of time-averaged densimetric Froude numbers $\langle v_F/ U_H \rangle _{\tau _{over}}$, where cases are categorized by the non-dimensional particle diameters $d_p/H$. For the exact values, see table 4. These show experiments to generally have a larger front speed than the simulations with this difference being more prominent when $d_p/H$ is small.

Table 4.

Experimental and simulation values of front velocities and averaged buoyant heights. The averages are taken over $\tau _{over}$, which is the time interval for which we have both experimental and simulation data in figures 10 and 12, respectively.

$^{a}$Simulation S0L15A was carried out with no-slip boundaries at $z=0, L_z$.

$^{b}$Simulations S0L15C and S3L15C were performed with $Sc=1$, as compared to the standard value of $Sc=7$.

The experimental data again show the pronounced seiche oscillations, while the simulation results exhibit somewhat smaller, more random-like variations. Nevertheless, experimental and numerical front velocities are seen to fluctuate around comparable mean values, and they show similar rates of decline. For a quantitative comparison, we average the experimental and computational front velocities over the time interval $\tau _{over}$ from $2 < t/T_H < 28$, in order to obtain $\langle v_F/ U_H \rangle _{\tau _{over}}$ shown in figure 10(b). The experimental front velocities are generally slightly higher than the computational values. This difference tends to increase as $d_p/H$ decreases, such that the difference is largest for the smooth case. We need to keep in mind, however, that for the smooth case the experimental Reynolds number was about twice that of the corresponding simulation, which may account for much of this difference.

In order to assess the influence of $Re$ and $Sc$ on the front velocity, we carried out two additional simulations for the smooth case. While S0L15B uses the standard set-up, S0L15A employs no-slip sidewalls, and S0L15C has a reduced Schmidt number value of $Sc=1$. We found that the lower $Sc$-value slightly reduced the front velocity, while the no-slip sidewall affected the front velocity by much less than one per cent.

3.4. Current height

The velocity of a gravity current, along with its destructive potential, are largely functions of the buoyant height. This height, in turn, depends on the net exchange of fluid between the frontal region and the tail section of the current, as well as on the leakage of current fluid into the porous bed. In this section, we will present results for the buoyant height of the current front as a function of time, and subsequently analyse the respective fluid exchange of the front with the tail section and the bed.

Let us consider the average height $h_M$ of a current within a distance $4H$ behind the front. Figure 11 presents experimental results for $h_M/H$ as a function of time for different parameter combinations. While all experiments show a decline of $h_M/H$ with time, the rate of this decline depends on the specific experimental conditions. Panel (a) demonstrates that, when the number of sphere layers is held constant at three, the buoyant height decreases more slowly in deeper water, i.e. for a larger ratio of buoyant height to bed height, $h_M/\eta$. In (b), we recognize that, in deep water, one, two and three layers of spheres result in approximately identical buoyant heights with similar rates of decline. However, this rate of decline is significantly larger than for a current over a flat bottom (E0L20, orange line).

Figure 11.

The spatially averaged current height $h_M$ as a function of time for four groups of experiments: (a) constant number of particle layers, but varying water depth; (b) constant water depth, but different numbers of layers; (c) water depth of 15 cm, varying number of layers; and (d) water depth of 10 cm, varying number of layers.

Panels (c,d) suggest that, for shallower total water depth, the difference between currents propagating over one, two or three layers increases. Taken together, these observations indicate that currents over porous beds behave in a fundamentally different way from those over flat beds, and that the exact depth of the porous bed becomes influential only in relatively shallow water. This suggests that a thick current traveling in deep water moves so fast that only fluid from the uppermost pores gets mixed into the current front. Shallower currents that travel more slowly, on the other hand, are more likely to allow fluid from the lower pores to rise and become mixed into the current front. The fluid exchange between the current and the bed will be analysed in more detail below.

Figure 12 compares the buoyant front height $h_M$ for the simulations to the corresponding experimental values. The time-averaged values $\langle h_M \rangle _{\tau _{over}}$ are provided in table 4. We note that, early on, the simulations tend to have larger buoyant front heights than the corresponding experiments, whereas the agreement improves at later times. This may at least partially be due to the perturbations introduced in the experiments by the removal of the gate and by the step in the sediment bed height, which are likely to promote an earlier transition to a fully turbulent state. Furthermore, the current height is seen to decrease for larger Reynolds number in the shallow simulations, which may again reflect the earlier transition to turbulence, cf. simulations S3L10A and S3L10B in figure 12(c). Finally, (d) shows that for smooth currents a lower Schmidt number is associated with a decrease in $h_M/H$, while for rough beds panel (b) shows no significant difference between the simulations for $Sc=7$ and $Sc=1$.

Figure 12.

Comparison of experimental and numerical results for the front height $h_M$. The agreement generally improves with time, which suggests that the discrepancy is mainly due to differences in the initialization of the currents.

The decrease in $h_M$ indicates that the current head loses salinity either into the bed or into the tail section of the current. This issue will be discussed in more detail below. In order to further analyse the dilution of the current front, we calculate the height $h_r$ as the maximum $y$-value where $\langle c\rangle _z=0.1$, at the fixed streamwise location $x=-2H$ relative to the front location of the current. Figure 13(a) shows that for the three-layer experiments and the smooth bottom case, respectively, $h_r$ decreases with time. Figure 13(b) indicates that $h_M$ decreases even more rapidly than $h_r$, so that the currents not only see a reduced height, but also become more dilute.

Figure 13.

Experimental measurements for $h_r/H$ (a), and $h_M/h_r$ (b) indicate that as their height decreases, the currents also become more dilute.

3.5. Streamwise velocity profiles

Figure 14 compares the streamwise velocity profiles for experiments and simulations of gravity currents over a smooth wall, as well as over one or three layers of spherical particles. The velocity profiles are evaluated at a distance $2H$ behind the front, and they are averaged both in the spanwise direction and over the time interval $13 \le t/T_H \le 15$. The peak velocities in both simulations and experiments are somewhat larger for the smooth wall, suggesting that rough sediment beds enhance the vertical mixing of streamwise momentum. Near the smooth wall, the velocity is smaller than near the surface of the sediment beds, which indicates that gravity currents propagating over rough beds see a small effective slip velocity at the bed surface. The above velocity profiles are similar to those observed by Zhou et al. (Reference Zhou, Cenedese, Williams, Ball, Venayagamoorthy and Nokes2017) in what the authors refer to as the overflowing regime.

Figure 14.

Streamwise velocity profiles, averaged in the spanwise direction, for zero, one and three layers of spheres. The profiles are evaluated at location $x= x_F-2H$, and they are averaged over the time interval $13 \le t/T_H \le 15$. When employing the PTV technique, the velocity values within approximately $3$ mm of the bottom and the top of the flume are affected by reflections, so that we do not show any experimental data in those regions.

4.1. Mass transfer

The experimental and simulation results presented in the preceding section show several interesting findings regarding the influence of the porous particle bed on the overall gravity current properties. Among them are the exchange of fluid between the current and the bed, a somewhat reduced front velocity that depends mostly on the particle size, a more rapid decrease of the current height and front velocity with time, and a modified velocity profile within the current due to enhanced vertical mixing of streamwise momentum. In order to obtain more detailed insight into the mechanisms by which the particle bed effects these changes, we now proceed to analyse the mass and momentum transfer between the current and the bed in depth.

We begin by analysing the transport of salinity into and out of the current head, in the reference frame moving with the front. Towards this end, we focus on the control volume $M(t)=[x_F(t)-x_{ML}, x_F(t)+ x_{MR}] \times [\eta , L_y] \times [0, L_z]$, with $x_{ML}=4H$ and $x_{MR}=0.5H$. Integrating transport equation (2.12) over this volume and dividing by its size $V_M= M_x H L_z$, where $M_x=x_{ML}+x_{MR}= 4.5H$, yields

(4.1)\begin{equation} \partial_t \frac{1}{V_M} \int_M c \, \textrm{d}V = \frac{1}{V_M} \int_{\partial M} \left [ v_F c \boldsymbol{e}_x - \boldsymbol{u} c + Pe_L^{-1} \boldsymbol{\nabla} c \right ]\boldsymbol{\cdot} \boldsymbol{n} \, \textrm{d}A ,\end{equation}

with $\boldsymbol {n}$ representing the outward normal to the boundary $\partial M$ of control volume $M$. By introducing the streamwise velocity in the moving reference frame $u_x^+= u_x-v_F$ and neglecting the small contributions due to diffusion, we obtain

(4.2)\begin{equation} \partial_t \langle c \rangle_M = \langle u_x^+ c^+ \rangle_{\partial M, W }\frac{1}{M_x} + \langle u_y^+ c^+ \rangle_{\partial M, S }\frac{1}{H}.\end{equation}

Here, the volume average over $M$ is denoted by $\langle \rangle _M$, while $\langle \rangle _{ \partial M, W}$ and $\langle \rangle _{ \partial M, S}$ indicate the surface averages of the fluxes across the western and southern boundaries, used to define the boundaries at $x^+ = -4H$ and $y=\eta$, respectively.

Figure 15 presents the temporal variation of the terms in the mass balance equation (4.2), for one flow over a smooth wall and three flows over three particle layers of different size. We observe that over time all current heads lose salinity, although the loss rate tends to decrease with time for all flows except the one with the largest particles. The relatively large temporal fluctuations in the rate of loss to the current tail may be related to vortical structures as well as variations in the interface height. Interestingly, the loss of current fluid into the bed increases with particle size $d_p/H$, whereas the loss of salinity to the tail (blue dashed-dotted line) decreases. The reasons behind this influence of the particle size will be discussed below.

Figure 15.

Mass balance terms according to (4.2), averaged over a moving window of two time units. The loss of current fluid into the bed (green line) increases with particle size, while the loss of head fluid to the tail (blue dash-dotted line) decreases for larger particles. $(a)\ \textrm{S0L15B}, d_{p}/ H = 0; (b)\ \textrm{S3L20A}, d_{p}/H = 0.056; (c)\ \textrm{S3L15B}, d_{p}/H = 0.081; (d)\ \textrm{S3L10B}, d_{p}/H = 0.14.$

Figure 16 presents simulation results in the moving reference frame for the averaged salinity flux $\langle c^+ u_x^+\rangle _{z,\tau }$ across the western boundary at $x^+ = -4H$. The profiles shown are averaged over a time interval $\tau =[\tau _0, \tau _1]$ that excludes the initial laminar evolution, with the specific time values provided in table 3. Compared to the flow over a smooth wall, the flows with particle layers exhibit lower salinity fluxes from the head into the tail just above $y = \eta$, and from the tail into the head further away from the lower boundary. This observation is consistent with the velocity profiles shown in figure 14, and it reflects the influence of the effective slip velocity at $y = \eta$ mentioned earlier for the flows over particle beds, along with increased vertical mixing in those flows. We will further explore the vertical mixing of momentum in § 4.2.

Figure 16.

Streamwise salinity flux in the moving reference frame $\langle u_x^+c^+ \rangle _{z,\tau }$ across the western boundary at $x^+=-4H$ as a function of the vertical coordinate, shown on a log scale. Compared to a flow over a smooth wall, the flows with particle beds exhibit reduced salinity fluxes from the head into the tail near the lower boundary, and from the tail into the head further away from the boundary.

Figure 17 analyses the time-averaged salinity flux from the current head into the bed, again in the moving reference frame. At the level of the bed surface, $y = \eta$, (a) shows that this flux has a pronounced maximum over the interval $0>x^+/H>-1$, while it levels off thereafter. This indicates that the current front is very effective at sweeping the ambient fluid out of the primary pore spaces and replacing it with current fluid. Here we define the primary pore spaces as the exposed pore spaces between the upper halves of the topmost layer of spheres. Larger particle sizes are associated with a larger primary pore space volume, and thus result in an increased vertical flux. The high-frequency oscillations reflect the periodicity of the bed arrangement.

Figure 17.

Time-averaged vertical salinity flux from the current head into the bed. (a) Flux at the top of the bed, $y = \eta$, and (b) flux at the centre of the first particle layer, $y= \eta - d_p/2$. The different lines correspond to simulations S3L20A, S3L15B, S3L10B, S3L10C and S1L15A (top to bottom of the legend). The brown short–long dashed line corresponds to a single layer of particles, whereas all other simulations have three layers. The small salinity flux into the bed for $x>0$ results from our definition of the front as the location where the current height $h=0.1$, so that the current extends slightly beyond $x=0$.

The subsequent transport of current fluid into the deeper pore layers occurs more slowly. This is demonstrated by figure 17(b), which shows that the vertical flux at the centre of the first particle layer, $y= \eta - d_p/2$, ramps up only gradually behind the current front. A comparison of (a,b) indicate that for $x/H < -1$ this flux into the deeper pore spaces roughly equals the flux at the bed surface, which suggests that by then the ambient fluid has largely been removed from the primary pore volume. The effective permeability of the bed increases with the particle size, which explains the enhanced salinity flux into the deeper layers of the bed for larger particles, as seen in (b). Interestingly, the $Re_H = 11\ 000$ simulation (beige solid line) transports more salinity into the bed than the $Re_H = 8000$ one, which may reflect the mobility increase due to the lower viscosity resistance. The single-layer simulation (brown line) shows a reduced flux due to a lower volume of bed fluid that can be exchanged, as compared to an equivalent three-layer simulation (green line).

From geometric considerations, it follows that the primary pore space volume $V_{inter}$ per $x,z$-base area $A$ is

(4.3)\begin{equation} V_{inter}/A=K_{inter}d_p, \quad\text{where } K_{inter}= \frac{\sqrt{3} - {\rm \pi}/3}{2\sqrt{3}}=0.1977 , \end{equation}

so that the current front overruns primary pore space volume per unit width at the rate $\langle v_F \rangle _{\tau } K_{inter}d_p$. Within the distance $H$ behind the front, current fluid enters the bed at the rate $\int _{-H}^{0} \langle u_y^+ c^+ \rangle _{z, \tau } \textrm {d}x$. Hence

(4.4)\begin{equation} \sigma= \frac{\int_{-H}^{0} \langle u_y^+ c^+ \rangle_{z, \tau} \textrm{d}x }{- \langle v_F \rangle_{\tau} K_{inter}d_p},\end{equation}

denotes the fraction of primary pore space fluid that is swept out within the distance $H$ behind the front. Figure 18 indicates that $\sigma$ increases with the particle diameter, and typically lies in the range of $0.5\pm 0.1$. This confirms our earlier hypothesis that a significant fraction of the primary pore fluid gets flushed out within a distance $H$ behind the front.

Figure 18.

Fraction $\sigma$ of the primary pore space fluid that is flushed out by the current within a distance $H$ behind the front.

4.2. Momentum transfer

To analyse the momentum budget, we integrate the momentum balance

(4.5)\begin{equation} \partial_t \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{\cdot} ( \boldsymbol{u} \boldsymbol{u})= -\boldsymbol{\nabla} p + Re^{-1} \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^T \right ) - \boldsymbol{e}_y c \end{equation}

over the moving control volume $\hat M=[x_F-4H, x_F+ 0.5H ]\times [\eta , \eta +\zeta ]\times [ 0, L_z]$ and divide by the spanwise domain width $L_z$. Note that $\hat M$ does not cover the entire vertical extent of the domain, but only the region from the bed surface at $y=\eta$ to the distance $\zeta$ above the bed where the streamwise velocity in the moving reference frame vanishes. Hence $y= \eta + \zeta$ can be viewed as the boundary between the current and the ambient counterflow. Figure 19 displays the spanwise-averaged value of $\zeta$ at different times, relative to the front position. It shows that both for a smooth wall and a rough bed $\zeta /H$ fluctuates around the value of 0.3 in the tail, so that we pick $\zeta /H=0.3$ as the upper boundary of the control volume for the momentum balance.

Figure 19.

Boundary between the current and the counterflow for evaluating the momentum balance, as defined in the text. (a) Flow simulation S0L15B over a smooth wall, and (b) flow simulation S3L20A over a bed of three particle layers.

In compact notation, we can write the $x$-component of the integrated (4.5) as

(4.6)\begin{equation} \zeta M_x \partial_t \langle u_x \rangle_{ \hat M} =\zeta { \omega}_f + \zeta {\tau}_p + M_x {\tau}_h+ M_x{ \tau}_b , \end{equation}

where the individual terms on the right-hand side are defined as

(4.7)\begin{equation} \left.\begin{array}{c@{}} \omega_f(t) = \langle u_x u_x - v_F u_x \rangle_{\hat M, W} - \langle u_x u_x - v_F u_x \rangle_{\hat M, E} , \\ \tau_p(t)= \langle\, p \rangle_{\hat M, W} - \langle\, p\rangle_{\hat M, E} , \\ \tau_h(t)= -\langle u_y u_x \rangle_{\hat M, N} + \dfrac{1}{Re_L} \langle \partial_y u_x \rangle_{\hat M, N} , \\ \tau_b(t)= +\langle u_y u_x \rangle_{\hat M, S}- \dfrac{1}{Re_L} \langle \partial_y u_x \rangle_{\hat M, S} . \end{array}\right\} \end{equation}

The terms in (4.6) represent, from left to right, the time rate of change of $x$-momentum within the control volume, the net convective momentum flux across the east and west boundaries ($\omega _f$), the net pressure stress ($\tau _p$), the momentum flux through the top boundary of the current at a distance $\zeta$ above the bed ($\tau _h$) and the bed stress ($\tau _b$) due to viscous stress and the convection of momentum. Note that we neglect the viscous stresses at the west and east boundaries, along with the $x$-derivative of the vertical velocity at the south and north boundaries, as these terms are small.

Figure 20 displays the terms of (4.6) for the representative simulation S3L10B. We observe that the bottom friction $\tau _b$ represents the main stress retarding the current, with a smaller contribution coming from the top friction $\tau _h$. These retarding forces are largely balanced by the convective momentum inflow $\omega _f$ and the streamwise pressure drop $\tau _p$, although the current still experiences a slight overall net deceleration. A detailed comparison with results from other simulations (not shown) indicates that the relative magnitude of the terms is similar for all simulations, although the absolute values can vary significantly, as will now be discussed for the bottom friction term.

Figure 20.

Contributions to the streamwise momentum balance in the control volume for simulation S3L10B, according to (4.6). The data have been smoothed by employing a moving average extending over two time units. The retarding influence of the bottom and top friction is largely balanced by the convective inflow of momentum and the streamwise pressure drop.

4.2.1. Bottom friction

In order to explore the dependence of the bottom friction on the particle size, we define the friction coefficient $C_f$ as

(4.8)\begin{equation} C_f(x,z,t) = -\frac{u_y u_x - Re_L^{-1} \partial_y u_x}{v_F^2/2} \quad \text{for } y=\eta . \end{equation}

Here, we included the convective transport of momentum as well, because a potential simplified model of the bed might assume an impermeable surface at $y=\eta$ and then $C_f$ could be used to represent the full momentum transport.

Time-averaged friction coefficient data (the averaging interval is given in table 3) are provided in figure 21. These demonstrate that the friction coefficient increases with the particle size. Even for the smallest particles, the bottom friction is substantially larger than for the smooth wall. Closer inspection of the results shows that most of this increase stems from the convective term $\langle u_y u_x \rangle _{z,\tau }$, which increases significantly for the largest particles with $d_p=0.1$.

Figure 21.

Time-averaged bottom friction coefficient $\langle C_f \rangle _{z, \tau }$ according to definition (4.8), for simulations S0L15B, S3L20A, S3L15B, S3L10C and S1L15A (from top to bottom in the legend).

Figure 22(a) displays the mean vertical velocity at the top of the bed as a function of the streamwise coordinate for selected simulations. The figure demonstrates that ahead of the front fluid is ejected out of the bed, while it is pushed into the bed below the current head. Panel (b) shows the variance of the vertical velocity at $y=\eta$ for the same set of simulations sharing the line style with (a). We find that its value increases with the particle size and the Reynolds number, as a result of the increased fluid mobility within the bed.

Figure 22.

(a) Mean vertical velocity component at the top of the bed for five simulations of flows over particle beds (S3L20A, S3L15B, S1L15A, S3L10B, S3L10C) (from top to bottom in the legend). Fluid is ejected out of the bed ahead of the front, and into the bed below the gravity current head. (b) Variance of the vertical velocity. This quantity increases with the particle size and the Reynolds number, reflecting the increased mobility of the fluid within the bed.

Based on Darcy's law, we can estimate the magnitude of the vertical velocity $u_D$ into the porous bed as

(4.9)\begin{equation} u_D \frac{\mu}{\tilde K}= \boldsymbol{\nabla} p, \end{equation}

where $\tilde K$ denotes the permeability of the bed. Carman (Reference Carman1997) showed empirically for flow through glass spheres with a porosity $\phi _f$ that the permeability can be approximated by

(4.10)\begin{equation} \tilde K = \frac{d_p^2\phi_f^3}{180(1-\phi_f)^2}. \end{equation}

For a densely packed bed with $\phi _f = 0.26$ and $d_p=0.01$ m, this results in a permeability of $\tilde K = 1.8\times 10^{-8}$ m$^2$. For a current height of approximately half the depth and negligible dilution, we can write (4.9) as

(4.11)\begin{equation} u_D =\frac{\frac{1}{2}{\rm \Delta} \rho g \tilde{H} \tilde{K}}{\rho_1\tilde{H}\nu} = \frac{1}{2}\mathcal{U}_H\frac{\tilde K}{\tilde H^2}Re_H,\end{equation}

or

(4.12)\begin{equation} u_D/ \mathcal{U}_H = \frac{1}{2}\frac{\tilde K}{\tilde H^2}Re_H. \end{equation}

We find that this predicts the correct order of magnitude for the downward velocities in figure 22(a). However, it significantly underpredicts the increasing velocities into the bed with increasing $d_p$, and overpredicts the increase in velocity with Reynolds number. The imperfections in applying the Darcy model to the flow into the bed are not unexpected because of the shallow depth of the rough bed, the relatively large size of the grains, and the sidewalls adding complexity to the flow. Due to the large pore size a Darcy–Brinkman formulation may be more applicable to the problem. The empirical constants in such a model are less accurately predictable and the model would still experience imperfections due to the sidewalls and the shallow depth of the bed.

Based on classical investigations of turbulent flows over rough walls (Nikuradse Reference Nikuradse1931; Schlichting Reference Schlichting1936; Pope Reference Pope2000), we expect the bottom friction to depend on the ratio of the particle radius to the thickness of the viscous sublayer. In order to be able to estimate the latter, we introduce the friction velocity

(4.13)\begin{equation} u_\tau =\frac{\sqrt{ \nu \partial_{ \tilde y} \tilde u_x}} {\mathcal{U}_L }= \sqrt{\partial_y u_x/ Re_L } , \end{equation}

and a length scale $L_W$ for the near-wall region

(4.14)\begin{equation} L_W = \frac{\nu}{\tilde u_\tau \tilde L_y} = (\partial_y u_x Re_L )^{-1/2} = u_\tau^{-1} Re_L^{-1}. \end{equation}

Figure 23(a) displays streamwise velocity profiles in wall units (solid line) at $x-x_F= -2H$, for three different simulations. A comparison with the linear profiles $u_\tau y/L_W+ u_x(y=\eta ) / u_\tau$ (dashed line) demonstrates that the viscous sublayer extends to the typical value of five wall units, i.e. $y-\eta \approx 5L_w$. Panel (b) shows the mean friction coefficient as a function of the particle radius normalized by the wall unit $L_W$. All of the cases with particle layers are seen to fall into the transitional roughness regime $10<d_p/(2L_W)<30$, indicating that the particle radius is slightly to substantially larger than the viscous sublayer thickness. Interestingly, the friction coefficient strongly increases once the particle radius exceeds about twenty wall units, as a result of the increased advective momentum mixing seen earlier in figure 22. This observation is consistent with the work of Fang et al. (Reference Fang, Han, He and Dey2018) and Manes et al. (Reference Manes, Pokrajac, McEwan and Nikora2009), who investigated open channel pressure-driven flows over fixed beds of spheres, finding enhanced momentum mixing for larger bed permeabilities.

Figure 23.

(a) Spanwise- and time-averaged streamwise velocity profiles, evaluated $2H$ behind the front and scaled with wall units, cf. (4.13) and (4.14). Panel (b) shows the averaged friction coefficient, (4.8), as a function of the particle size scaled with the wall length $L_W$. The values of the wall lengths are $L_W=2.2\times 10^{-3}$, $2.2\times 10^{-3}$, $1.9 \times 10^{-3}$, $2.2\times 10^{-3}$, $1.8\times 10^{-3}$ and $2.0\times 10^{-3}$, respectively, for simulations S0L15B, S3L20A, S3L15B, S3L10B, S3L10C and S1L15A. Note that the friction coefficient increases significantly when the particle radius exceeds approximately twenty wall units.