2. Description of the numerical model and simulations

The filtered dimensionless Navier–Stokes equations are solved along with the advection–diffusion equation for the non-dimensional scalar concentration on a non-uniform Cartesian mesh (Ooi, Constantinescu & Weber Reference Ooi, Constantinescu and Weber2007). The dimensionless filtered concentration is $C=(C^{\prime }-C_{min}^{\prime })/(C_{max}^{\prime }-C_{min}^{\prime })$ , where $C^{\prime }$ is the concentration, and $C_{max}^{\prime }$ and $C_{min}^{\prime }$ represent the maximum and minimum initial concentrations in the domain. The concentrations $C$ and $C^{\prime }$ are linearly related to the density of the fluid. The reduced gravity is denoted by $g^{\prime }=g({\it\rho}_{max}-{\it\rho}_{min})/{\it\rho}_{max}$ , where ${\it\rho}_{max}\equiv {\it\rho}+{\rm\Delta}{\it\rho}$ and ${\it\rho}_{min}\equiv {\it\rho}$ are the initial maximum (lock fluid) and minimum (ambient fluid) densities in the domain, and $g$ denotes gravitational acceleration. The channel height, $H$ , and the buoyancy velocity, $u_{b}=\sqrt{g^{\prime }H}$ , are used as the characteristic length and velocity scales respectively, with which the variables are non-dimensionalized and times are rendered dimensionless in terms of $t_{0}=H/u_{b}$ . It is assumed that the density differences are sufficiently small ( ${\rm\Delta}{\it\rho}/{\it\rho}\ll 1$ ) that the Boussinesq approximation may be invoked and they are only significant in the gravitational terms of the momentum equations. The subgrid-scale viscosity, ${\it\nu}_{t}$ , and subgrid-scale diffusivity, ${\it\kappa}_{t}$ , in the filtered momentum equations and in the advection–diffusion equation for $C$ are calculated using a dynamic Smagorinsky model. The equations used to evaluate ${\it\nu}_{t}$ and ${\it\kappa}_{t}$ are given in Pierce & Moin (Reference Pierce and Moin2001). The use of a dynamic model eliminates the need to specify any stratification correction in the subgrid-scale model (Kirkpatrick et al. Reference Kirkpatrick, Ackerman, Stevens and Mansour2006; Rodi, Constantinescu & Stoesser Reference Rodi, Constantinescu and Stoesser2013). The dynamic model was found to be very successful in predicting scalar transport and mixing in turbulent flows, especially for cases where the resolved scales drive the mixing process (Rodi et al. Reference Rodi, Constantinescu and Stoesser2013). The two parameters in the non-dimensional governing equations are the channel Reynolds number, $\mathit{Re}=u_{b}H/{\it\nu}$ , and the molecular Schmidt number, $\mathit{Sc}={\it\nu}/{\it\kappa}$ , where ${\it\nu}$ is the kinematic viscosity of the fluid and ${\it\kappa}$ is the molecular diffusivity of the dissolved species that gives rise to the density field. The coordinates in the three directions are denoted either ( $x_{1},x_{2},x_{3}$ ) in index notation or ( $x,y,z$ ). The vertical axis is $y$ , aligned with gravitational acceleration, while the $x$ -axis is aligned along the channel axis and the $z$ -axis is spanwise.

The finite-volume LES code (Pierce & Moin Reference Pierce and Moin2001; Chang & Constantinescu Reference Chang and Constantinescu2013) advances the governing equations in time using a semi-implicit iterative method that employs a staggered conservative space–time second-order-accurate discretization. A Poisson equation is solved for the pressure using a multigrid approach. All operators are discretized using central discretizations, except the convective term in the advection–diffusion equation for the concentration of the dissolved species, for which the quadratic upwind interpolation for convective kinematics (QUICK) scheme (Leonard Reference Leonard1979) is used. Despite the non-monotonic properties of the scheme that allow potentially non-physical values of the scalar to be generated in the simulation, the present code was found to predict accurately scalar transport in turbulent flow applications involving both passive and non-passive scalars (see, for example, the detailed validation for combustion applications discussed by Pierce & Moin (Reference Pierce and Moin2001, Reference Pierce and Moin2004)). Although the QUICK scheme is dissipative, for high-Reynolds-number flows the dissipation introduced by the subgrid-scale model is significant. This is the main reason why the numerical diffusion associated with the QUICK scheme in LES of scalar transport does not have a very adverse effect on the quality of the solution. Chang, Constantinescu & Park (Reference Chang, Constantinescu and Park2006, Reference Chang, Constantinescu and Park2007) discussed validation of the numerical method to predict ejection of non-buoyant and buoyant contaminants from bottom channel cavities. Ooi et al. (Reference Ooi, Constantinescu and Weber2007, Reference Ooi, Constantinescu and Weber2009) performed LES of lock-exchange intrusion currents and gravity currents propagating over a smooth flat surface. Comparison with laboratory experimental realizations showed that 3D LES accurately captured the structure of the current at different stages of its evolution. Large eddy simulation correctly predicted the change in the front speed, $U_{f}$ , with Reynolds number, Re, during the slumping phase for values up to $10^{6}$ and that $U_{f}\sim t^{-1/3}$ during the buoyancy-inertia phase, in accord with analytical models of the motion (Hoult Reference Hoult1972). Gonzalez-Juez, Meiburg & Constantinescu (Reference Gonzalez-Juez, Meiburg and Constantinescu2009) and Gonzalez-Juez et al. (Reference Gonzalez-Juez, Meiburg, Tokyay and Constantinescu2010) used the same numerical model to study the changes in structure of the gravity current as it interacts with isolated obstacles mounted on a flat bed surface or situated at a small distance from the channel bed. Large eddy simulation successfully reproduced the measured time-varying drag and lift coefficients of gravity current flows impinging circular and rectangular cylinders (Ermanyuk & Gavrilov Reference Ermanyuk and Gavrilov2005a ,Reference Ermanyuk and Gavrilov b ). Tokyay, Constantinescu & Meiburg (Reference Tokyay, Constantinescu and Meiburg2011, Reference Tokyay, Constantinescu and Meiburg2012, Reference Tokyay, Constantinescu and Meiburg2014) showed that LES predictions of vertical profiles of mean velocity and turbulent kinetic energy agreed well with experimental measurements conducted for turbulent lock-exchange gravity currents over a smooth flat surface and that the front speed variation during the drag-dominated regime was in close agreement with shallow flow theory for lock-exchange currents propagating over a rough surface.

All lock-exchange simulations reported in the current study were conducted in a channel with smooth horizontal top and bottom walls and lateral vertical endwalls. The channel contained an array of square cross-section cylinders of side length $D$ . Their location was close to a regular staggered pattern in the $x{-}y$ planes, but a small random displacement was applied to the initial regular pattern to reduce the regularity of the wake-to-wake interactions between neighbouring cylinders. The magnitude of the random displacement in the vertical and horizontal directions varied between $0.5D$ and $1.0D$ . No-slip boundary conditions were employed on the velocity field at the channel walls and on the surfaces of the cylinders. The surface-normal concentration gradient was set to zero at all wall boundaries. The flow was assumed to be periodic in the spanwise direction. The time step was $0.001t_{0}$ . The Schmidt number was 600, which corresponds to saline diffusion in water. Härtel et al. (Reference Härtel, Meiburg and Necker2000), Necker et al. (Reference Necker, Härtel, Kleiser and Meiburg2005) and Cantero et al. (Reference Cantero, Lee, Balachandar and Garcia2007) showed that variations in the value of the Schmidt number did not strongly alter the results in lock-exchange simulations of gravity currents as long as it was of order unity or larger. This conclusion was confirmed by Ooi et al. (Reference Ooi, Constantinescu and Weber2009) using the same code as that employed to perform the present simulations.

The numerical runs were performed with channel Reynolds numbers, Re, in the range 100–150 000. Since a typical temperature difference between shaded and unshaded regions is just a few degrees Celsius, in a typical water depth of a few metres, it is implied that exchange flows within aquatic canopies are characterized by a Reynolds number of order $10^{4}$ . Furthermore, this Reynolds number is typical of many laboratory realizations of these flows. However, in other applications in which arrays of porous barriers are used to decelerate currents, the Reynolds number may be rather higher (e.g. $\mathit{Re}=10^{5}{-}10^{8}$ ), because the typical length scales and flow speeds are increased (see Naaim-Bouvet et al. Reference Naaim-Bouvet, Naaim, Bacher and Heiligenstein2002; Oehy & Schleiss Reference Oehy and Schleiss2007). However, as will be shown below, the motion is retarded by the obstacles and decelerates, thus promoting the relative magnitude of viscous processes. Eventually the motion becomes fully controlled by viscous effects, and for this reason we have also conducted simulations at relatively low Reynolds numbers so that the effects of fluid turbulence are suppressed and the viscous dynamics may be analysed.

The obstructed channel is geometrically characterized by the following variables. There are $N$ square cylinders within the domain of volume $V$ , so the volume fraction of solids is given by

where $W$ is the width of the domain. The average spacing between the cylinders, $s$ , is given by

Finally, the frontal area of the cylinders per unit volume, $a$ , is given by

The inverse of this quantity, $1/a$ , is a length scale characterizing the porous medium.

The simulations were conducted with volume fractions of solids ${\it\phi}$ between 0 and 25 % (see table 1). Although most vegetated canopies generally have ${\it\phi}<8\,\%$ , there are practical applications where higher solid volume fractions can occur (Jamali et al. Reference Jamali, Zhang and Nepf2008). For example, in the case of arrays of porous barriers used to decelerate turbidity currents advancing in reservoirs, ${\it\phi}$ generally varies between 5 and 20 % and even in the case of vegetated canopies ${\it\phi}$ values as high as 45 % have been reported in mangroves with mean trunk diameters in the range $D=4{-}9~\text{cm}$ (Furukawa, Wolanski & Mueller Reference Furukawa, Wolanski and Mueller1997). The range of ${\it\phi}$ in the experiments of Zhang & Nepf (Reference Zhang and Nepf2008) for currents propagating in an aquatic canopy was 0–35 %. The range of $D/H$ was 0.032–0.07 in the present simulations. The average value of the cylinder Reynolds number during the quadratic drag-dominated regime, $\mathit{Re}_{D}$ , was between 60 and 250 for the $\mathit{Re}=15\,000$ simulation and close to 1700 for the $\mathit{Re}=150\,000$ simulation (table 1). These values are within the range considered in the experiments of King et al. (Reference King, Tinoco and Cowen2012) conducted for steady constant-density flow in a channel containing a layer of aquatic vegetation.

Table 1.

The main parameters of the numerical simulations. As in the text, $a$ is the frontal area of the cylinders per unit volume, $\mathit{Re}_{D}$ is the cylinder Reynolds number defined with $D$ and the mean $U_{f}$ after the transition to the drag-dominated regime has started, $N$ is the total number of cylinders, $C_{D}$ is the mean total drag coefficient estimated using the instantaneous velocity fields, $C_{D}^{\prime }$ is the mean total drag coefficient estimated using $U_{f}(t)$ as the velocity scale and ${\it\Gamma}_{D}=C_{D}aH/(1-{\it\phi})=C_{D}({\it\phi}H/D)/(1-{\it\phi})$ is a dimensionless combined drag parameter for the region containing the cylinders.

Apart from the cylinder diameter relative to the channel depth, $D/H$ , the other non-dimensional length scale is $aH$ . All present simulations were conducted with $aH>0.1$ (table 1) and with relative spacing, $s/D$ , greater than 2. As described below, we used the simulations to estimate the drag coefficient for the array of cylinders. This is denoted $C_{D}$ and is reported in table 1, along with a combined drag parameter for the porous medium, ${\it\Gamma}_{D}=C_{D}aH/(1-{\it\phi})=C_{D}{\it\phi}(H/D)/(1-{\it\phi})$ , which was larger than 0.48 (table 1). This latter parameter represents the ratio of obstacle drag force to momentum transport and can be also thought of as a global measure of the resistance provided by the cylinders; it will be shown to play a key role in the analysis that follows (§ 3). As discussed by King et al. (Reference King, Tinoco and Cowen2012), flows with ${\it\Gamma}_{D}>0.1$ and $s/D>1$ correspond to the ‘dense flow regime’ in which the size of the wake eddies and the turbulent kinetic energy generated by the cylinders scale with $D$ , and vertical gradients in the flow are fairly negligible. The range of ${\it\Gamma}_{D}$ in our simulations was 0.48–8.52, which overlaps with the values estimated by Jamali et al. (Reference Jamali, Zhang and Nepf2008) for field applications. We note that most experiments were conducted with ${\it\Gamma}_{D}<5$ , although ${\it\Gamma}_{D}$ values as high as 70 were recorded in their experiments with ${\it\phi}=35\,\%$ .

The computational mesh containing square cylinders was stretched in the wall-normal direction of no-slip surfaces (top and bottom walls, lateral endwalls) to resolve the near-wall flow. The grid resolution was also increased in the vicinity of the faces of the cylinders. In many of the simulations the grid spacing in these regions was approximately $0.0035H$ . Preliminary simulations were conducted for constant-density uniform flow past isolated cylinders. The grid resolution close to the isolated cylinder was the same as the one used for the cylinders within the array in the lock-exchange simulations, and the isolated cylinder simulations were conducted at comparable cylinder Reynolds numbers. The predictions of the drag coefficient for the isolated cylinder were found to be in good agreement with data in the literature (Yoon, Yang & Choi Reference Yoon, Yang and Choi2010). In most of the lock-exchange simulations, the grid contained $2400\times 200\times 48$ nodes in the streamwise (domain length $12H$ ), vertical (domain width $H$ ) and spanwise (domain height $H$ ) directions respectively. In the simulation conducted with $\mathit{Re}=150\,000$ (VHR_SVF12_D50), the grid was finer ( $4800\times 400\times 64$ ) because of the need to sufficiently resolve the near-wall flow. In the simulation with no cylinders (HR_SVF00) in which the current advances much faster, the number of grid points in the streamwise direction was 9600 to account for the increase in the channel length to $48H$ .

Despite the fact that the arrangement of the cylinders was not exactly reflectionally symmetrical, the temporal evolution of the main variables characterizing the evolution of the currents along the bottom and top of the domain was found to be close to anti-symmetrical. This is why only the evolution of the denser bottom-propagating current forming on the right side of the lock gate in figure 1 is analysed. Typical numerical results are shown in figures 2(c) and 3. These are discussed in detail below, but first we present a shallow layer model of the motion.

Figure 2.

(a) The height and (b) the depth-averaged velocity fields as a function of $x/x_{f}(t)$ for the similarity solution of a drag-dominated lock-exchange flow in the regime of high Reynolds number ( $C_{D}=\text{const.}$ ). Also plotted with symbols is the mean streamwise velocity determined from the numerical solution. (c) The scaled concentration field, $C$ , as a function of downstream distance, $x$ , and height, $z$ , when $x_{f}/H=4.6$ , from the simulation with $\mathit{Re}=15\,000$ , ${\it\phi}=12\,\%$ and $D/H=0.05$ . The dashed (blue online) and dash-dotted (red online) lines in panel (c) show the interface shape given by the analytical model, $h/H$ , and from the numerical simulation, $h_{n}/H$ , respectively. (d) The dimensionless position of the front, $x_{f}/H$ , as a function of dimensionless time $t/t_{0}$ , showing the similarity solution (dashed black line) and the results from the numerical simulation (solid pink line).

Figure 3.

Sketch of the front position, $x_{f}$ , as a function of time, $t$ , and the transitions between the various dynamical regimes for currents propagating into an obstructed channel predicted by shallow flow theory. The exponents are predicted by theory within different regimes (slumping phase, quadratic drag-dominated regime, linear drag-dominated regime). The transition times ( $t_{1},t_{2}$ ) and locations where transition is expected to occur ( $x_{1},x_{2}$ ) are defined in § 3. In frame (b) the current transitions directly from the slumping phase to the linear drag-dominated regime, and this corresponds to low-Reynolds-number lock-exchange flows.

3. Shallow water theory model

We analyse the density-driven motion to calculate the evolution of the interface between the dense and less dense fluids. We assume that mixing processes are negligible, so that the fluid density is ${\it\rho}+{\rm\Delta}{\it\rho}$ if it lies below the interface $(y=h)$ or ${\it\rho}$ if above, and the interface evolves according to a kinematic condition. In dimensional form, mass conservation in each layer is then given by

where $u_{1}$ and $u_{2}$ are the depth-averaged velocities of the fluid phase in the lower and upper layers respectively, and $y=0$ and $y=H$ are the lower and upper boundaries of the channel (see figure 1). We note that since there is no imposed external flux along the channel, there is volume flux balance given by

The motion is analysed in a regime where the velocity is predominantly horizontal; vertical fluid accelerations are negligible and so the pressure, $p$ , adopts a locally hydrostatic balance given by

where $p_{T}(x)$ denotes the as yet undetermined pressure along the upper boundary $y=H$ . Horizontal gradients of the pressure distribution drive the motion, while hydraulic resistance is provided by the array of cylinders through which the fluid flows. The drag at the boundaries of the porous domain is assumed to be much smaller than the drag due to the dispersed cylinders. The momentum balances for each layer in the Boussinesq regime are given by

(cf. Hatcher et al. Reference Hatcher, Hogg and Woods2000; Nepf Reference Nepf2012), where the drag is parameterized via the introduction of a drag coefficient $C_{D}$ . In a regime for which the cylinder Reynolds number is large ( $\mathit{Re}_{D}\gg 1$ ), the drag coefficient becomes independent of the Reynolds number, although it could still be a function of the other geometrical parameters specifying the porous domain ( ${\it\phi}$ and $D/H$ ). We note that if the Reynolds number is sufficiently small ( $\mathit{Re}_{D}\ll 1$ ), then $C_{D}\sim 1/\mathit{Re}_{D}$ and the hydraulic resistance becomes linearly proportional to the velocity. In this situation, the flow is analogous to motion through a porous medium modelled by Darcy’s law (Huppert & Woods Reference Huppert and Woods1995).

Following initiation from full-depth lock-release conditions, namely $h=H$ when $x<0$ and $h=0$ when $x>0$ , the dynamics governed by (3.5) and (3.6) is dominated by inertia, and the effects of drag induced by the cylinders have little effect on the motion. However, progressively the drag induced by the cylinders and the buoyancy force become the dominant terms in (3.5) and (3.6), while the inertial terms become negligible. It is insightful to employ scaling analysis to determine the temporal dependence of the velocity within the bottom layer, $u_{1}$ , and extent of the flow, $x_{f}$ , noting that the height of the interface $h\sim H$ for the lock-release flows considered here. Thus, when drag is negligible, we find $u_{1}\sim (g^{\prime }H)^{1/2}$ and $x_{f}\sim (g^{\prime }H)^{1/2}t$ . With these dependences, we may estimate the magnitude of the inertial terms in (3.5) as $(1-{\it\phi})\partial u_{1}/\partial t\sim (1-{\it\phi})(g^{\prime }H)^{1/2}/t$ , while the drag $C_{D}{\it\phi}u_{1}|u_{1}|/2D\sim C_{D}{\it\phi}g^{\prime }H/2D$ . They are thus comparable when $t\sim t_{1}=2(1-{\it\phi})D/[C_{D}{\it\phi}(g^{\prime }H)^{1/2}]$ and $x_{f}\sim x_{1}=2(1-{\it\phi})D/(C_{D}{\it\phi})$ . When $t\ll t_{1}$ , drag is negligible, but when $t\gg t_{1}$ , drag becomes a dominant force and changes the momentum balance. Importantly, we note that the dimensional dependences of these transition length and time scales may be written as

These expressions therefore reveal ${\it\Gamma}_{D}=C_{D}{\it\phi}D/[(1-{\it\phi})H]$ as the key dimensionless group of parameters that determines the onset of drag effects and their subsequent magnitude.

We now examine the regime in which buoyancy and drag forces balance, noting that the momentum equations for each layer simplify to balance between the streamwise pressure gradient and the hydraulic resistance due to the cylinder (see Hatcher et al. Reference Hatcher, Hogg and Woods2000):

We may then eliminate the drag forces using the flux balance (3.3) to deduce

Then we find that the momentum balance in the lower layer may be written as

This leads to the governing equation for the interface $h(x,t)$ ,

where $K=2g^{\prime }D(1-{\it\phi})/[C_{D}{\it\phi}]=2u_{b}^{2}/{\it\Gamma}_{D}$ , a constant of dimensions $L^{2}T^{-2}$ that characterizes the mobility of the density-driven flow. We seek a similarity solution to (3.12), by writing $h=HF({\it\eta})$ , with the similarity variable defined by

where $B$ is a dimensionless constant chosen such the interface reaches the lower boundary at ${\it\eta}=1$ ; this implies that $x_{f}(t)=B(KHt^{2})^{1/3}$ . Substitution of this form into the governing equation (3.12) leads to an ordinary differential equation (ODE) for $F({\it\eta})$ , given by

subject to boundary conditions $F(0)=1/2,F(1)=0$ and $F^{\prime }(1)=-4B^{3}/9$ . The first of these boundary conditions demands that the intrusions of dense and less dense fluid are symmetric; the second defines the front to occur at ${\it\eta}=1$ ; the third is forced by the singular nature of the governing equations at ${\it\eta}=1$ . This system is therefore a differential eigenvalue problem, which is solved straightforwardly using numerical techniques to evaluate the constant $B=1.162$ and to find the interface shape, $h/H$ (figure 2 a). Also plotted (figure 2 b) is the dimensionless velocity field in the lower layer, given by

Figure 2(c) shows the interface shape given by the analytical model superimposed on the concentration field predicted in a simulation with ${\it\phi}=12\,\%$ , $D/H=0.05$ and $\mathit{Re}=15\,000$ . Also shown is the interface shape, $h_{n}/H$ , calculated based on the concentration field predicted by the numerical solution, where $h_{n}=\int _{0}^{H}\!C\text{d}y$ . The overall agreement between the two curves in figure 2(c) is good, with the analytical solution predicting slightly lower values of the interface height, $h$ , over the downstream parts of the current. A similar level of agreement is observed in figure 2(b) between the numerical solution and the analytical model for the distribution of the non-dimensional mean streamwise velocity, $u_{1}/U_{f}$ . The agreement is satisfactory given that mixing is not accounted for in the analytical model. The agreement is worst close to the front of the current where mixing is especially strong. We also plot in figure 2(d) a comparison between the position of the front measured in the numerical simulations and the similarity solution ( $x_{f}/H=B(2/{\it\Gamma}_{D})^{1/3}(t/t_{0})^{2/3}$ ), noting that there is fair agreement between the two. In order to make this comparison and to estimate ${\it\Gamma}_{D}$ , we have determined the drag coefficient for the array from the simulation data using the mean velocity within the bottom current as the velocity scale (more details are given in § 6).

This solution is based upon the assumption that the drag exerted by the cylinders may be represented by a constant drag coefficient. Since the velocity is progressively slowing, this implies that there is a time scale, $t_{2}$ , at which the motion is no longer adequately modelled by the high-Reynolds-number representation of the drag from each cylinder. This occurs when $\mathit{Re}_{D}\sim 1$ , which is given when

The model of the motion developed above is appropriate then for $t_{1}\ll t\ll t_{2}$ , but for $t_{2}\ll t$ the governing equations are modified so that the drag is now linearly proportional to the flow speed. In this scenario the dominant balance between the hydrostatic pressure gradient and the drag is given by

where ${\it\gamma}$ is a dimensionless constant and the drag coefficient $C_{D}={\it\gamma}/\mathit{Re}_{D}$ . The governing system is analogous to flow within a porous medium governed by Darcy’s law. In this situation Huppert & Woods (Reference Huppert and Woods1995) establish similarity solutions for the motion given by

where $x_{f}(t)=(g^{\prime }H(1-{\it\phi})D^{2}t/[{\it\gamma}{\it\nu}{\it\phi}])^{1/2}$ .

To summarize, for these gravity currents initiated from full-depth lock-release conditions and propagating in an obstructed channel with a sufficiently high volume fraction of solids, the self-similar solution obtained based on shallow flow theory predicts (a) $U_{f}\sim t^{-1/3}$ and $x_{f}\sim t^{2/3}$ at sufficiently high Reynolds numbers for which one can assume $C_{D}$ to be constant, and (b) $U_{f}\sim t^{-1/2}$ and $x_{f}\sim t^{1/2}$ at sufficiently low Reynolds numbers where one can consider that a linear drag law ( $C_{D}\sim 1/\mathit{Re}_{D}$ ) applies (cf. Huppert & Woods Reference Huppert and Woods1995; Tanino et al. Reference Tanino, Nepf and Kulis2005). Furthermore, one expects currents to undergo transitions between behaviours. At very early times ( $t\ll t_{1}$ ), the motion is inertia dominated, before becoming drag dominated but with the drag proportional to the square of the velocity ( $t_{1}\ll t\ll t_{2}$ ), before entering a regime in which the drag is linearly proportional to the velocity ( $t_{2}\ll t$ ). This is illustrated in figure 3(a). It is possible that the intermediate regime of quadratic drag would not be observed (figure 3 b). This would occur if $t_{2}\ll t_{1}$ , which corresponds to the cylinder Reynolds number being small ( $\mathit{Re}_{D}\ll 1$ ).

4. Structure of the gravity current

5. Total streamwise drag force

The streamwise drag force on an individual cylinder is calculated by adding the streamwise components of the pressure and viscously derived forces, and then the total drag force on the cylinders $F_{D}$ is the sum of the magnitude of these forces throughout the entire computational domain. It should be noted that because the current evolves symmetrically (to a close approximation), $F_{D}/2$ is the magnitude of the force in each layer. The effects of the parameters ${\it\phi},D/H$ , $aH$ and Re on the total dimensionless streamwise drag force on the cylinders, $F_{D}/({\it\rho}u_{b}^{2}HW)$ , are discussed next.

The evolution of the streamwise drag $F_{D}$ for high-Reynolds-number currents (figure 9 a,b) is characterized by an initial regime in which $F_{D}$ increases sharply, followed by a transient regime which is characterized by a more gradual increase of $F_{D}$ and larger fluctuations about the mean value. The channel in the high-Reynolds-number simulations was not long enough to establish unequivocally that the flow reaches a regime in which $F_{D}$ fluctuates around a constant value, although the data are suggestive of this behaviour.

Figure 9.

The total drag force acting on the cylinders, $F_{D}/({\it\rho}u_{b}^{2}HW)$ . (a) The effect of the solid volume fraction in the $\mathit{Re}=15\,000$ simulations with ${\it\phi}=1.25\,\%$ , $D/H=0.032$ (light blue line L3), ${\it\phi}=5\,\%$ , $D/H=0.032$ (blue line L2) and ${\it\phi}=12\,\%$ , $D/H=0.032$ (dark grey line L1); (b) the effect of the relative cylinder size, $D/H$ , in the simulations with $\mathit{Re}=15\,000$ and ${\it\phi}=12\,\%$ and with $D/H=0.032$ (dark grey line L1), $D/H=0.05$ (pink line L4) and $D/H=0.07$ (dark blue line L5); (c) the effects of Re in the ${\it\phi}=12\,\%$ simulations with $\mathit{Re}=375$ , $D/H=0.05$ (orange line L9), $\mathit{Re}=15\,000$ , $D/H=0.05$ (pink line L4) and $\mathit{Re}=150\,000$ , $D/H=0.05$ (green line L6).

Simulations conducted with $\mathit{Re}=15\,000$ and $D/H=0.032$ show that as the solid fraction increases, the distance the front propagates before the start of the transient regime for $F_{D}$ becomes larger. At the same time, the magnitude of the fluctuations of $F_{D}$ increases with ${\it\phi}$ . This is because, as the distance between the cylinders increases, the coherence of the eddies detaching from the separated shear layers and of the eddies shed in the wake of the cylinders increases (figure 5 b–d). The increased coherence of such eddies increases the amplitude of the larger-scale fluctuations of the forces on the cylinder from which the eddy originated and also on the cylinders situated along its path. As all simulations in figure 9(a) were conducted with identically sized cylinders, the monotonic increase of $F_{D}$ with ${\it\phi}$ is due to the increase of the number of cylinders per unit area. Finally, the results in figure 9(a) for the simulation with ${\it\phi}=12\,\%$ suggest that the constant- $F_{D}$ regime is reached when $x_{f}/H\cong 4.5$ .

The effect of increasing the cylinder size, $D/H$ , or equivalently the frontal area, $aH={\it\phi}/(D/H)$ , while maintaining ${\it\phi}$ constant is to reduce $F_{D}$ , at least until the later stages of the transient regime when $F_{D}$ appears to approach the same constant value in the simulations conducted with $D/H=0.032$ and $D/H=0.05$ (figure 9 b). By contrast, the effect of increasing ${\it\phi}$ and $D/H$ while maintaining $aH$ constant does not result in significant changes in the levels of $F_{D}$ when plotted as a function of $x_{f}$ (e.g. compare the $\mathit{Re}=15\,000$ cases with ${\it\phi}=5\,\%$ , $D/H=0.032$ and with ${\it\phi}=12\,\%$ , $D/H=0.07$ in figure 9).

In the simulation conducted with $\mathit{Re}=375$ , the total drag force reaches a regime where $F_{D}$ is approximately constant after $t/t_{0}\cong 40$ when $x_{f}/H\cong 3$ (figure 9 c). The drag force also becomes approximately constant in the simulation with $\mathit{Re}=175$ (not shown). As no eddies separate behind the cylinders and the degree of unsteadiness of the separated shear layers is very small (figure 8 b), the amplitude of the large-scale fluctuations of $F_{D}$ in the simulation with $\mathit{Re}=375$ is very small compared with those observed in the simulations with $\mathit{Re}=15\,000$ and 150 000 conducted with the same values of ${\it\phi}$ (12 %) and $D/H$ (0.05). The increase of the Reynolds number from $\mathit{Re}=15\,000$ to $\mathit{Re}=150\,000$ does not result in significant differences in the mean variation of $F_{D}$ with $x_{f}/H$ .

The shallow water theory of § 3 predicts that $F_{D}/({\it\rho}u_{b}^{2}HW)$ should become constant once a balance between drag and buoyancy develops in the later stages of the motion. If $f_{D}$ is the drag force per unit volume due to the cylinders then

This equation holds for both high- and low-Reynolds-number currents in which $C_{D}=\text{constant}$ and $C_{D}\sim 1/\mathit{Re}_{D}$ respectively. We first proceed by assuming that $\mathit{Re}_{D}\gg 1$ and so

Conversely, when $\mathit{Re}_{D}\ll 1$ , we find that

In the regime of large Reynolds number, we find that this prediction for $F_{D}$ is closely borne out by the numerical simulations in which $F_{D}$ reaches a constant regime (e.g. the asymptotic value of $F_{D}/({\it\rho}u_{b}^{2}HW)$ was 0.24 in the $\mathit{Re}=15\,000$ simulation with ${\it\phi}=12\,\%$ and $D/H=0.032$ , which is close to the theoretical value of 0.25). The distance travelled by the current before the constant regime is reached is expected to increase with decreasing ${\it\phi}$ (see (3.7)). For example, the $\mathit{Re}=15\,000$ simulations with ${\it\phi}<10\,\%$ are still far from reaching the constant- $F_{D}$ regime by the time the current approaches the end of the channel $(x=6H)$ . In the regime of small Reynolds numbers, $\mathit{Re}_{D}\ll 1$ , we find that the measured force is rather smaller than predicted (e.g. $F_{D}/({\it\rho}u_{b}^{2}HW)$ predicted in the $\mathit{Re}=100$ and $\mathit{Re}=375$ simulations with ${\it\phi}=12\,\%$ is 0.22, while the theoretical value is 0.29).

6. Front velocity

In this section we analyse the unsteady propagation of the front position and its parametric dependence upon the volume fraction of cylinders, ${\it\phi}$ , the channel Reynolds number Re, and the combined drag coefficient for the porous medium, ${\it\Gamma}_{D}=C_{D}aH/(1-{\it\phi})$ . These were identified as the key parameters in the shallow layer modelling of § 3 and in Jamali et al. (Reference Jamali, Zhang and Nepf2008). However, before examining these dependences, we first use the $F_{D}$ estimated from the numerical simulations to deduce the drag coefficients of the cylinders.

7. Global energy budget

For compositional gravity currents, it is convenient to form the following differential equation to relate the rates of change of the potential and kinetic energies (Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2005; Ooi et al. Reference Ooi, Constantinescu and Weber2009):

In these expressions, the total kinetic energy and excess potential energy within the flow domain ${\it\Omega}$ are denoted by $E_{k}$ and $E_{p}$ respectively and ${\it\varepsilon}$ is the total dissipation rate. An integral balance equation for the mechanical energy can be obtained by integrating (7.1) with time,

where $E_{p0}$ is the total initial excess potential energy. The initial kinetic energy is equal to zero because the simulations are started from rest. The term $E_{d}=\int _{0}^{t}{\it\varepsilon}({\it\tau})\text{d}{\it\tau}$ is the time integral of the total dissipation rate. In all simulations, mechanical energy over the computational domain, expressed through (7.2), was conserved with an error of less than 0.5 % at all times.

In the simulations for obstructed channels (with ${\it\phi}>3\,\%$ ), after transition to the drag-dominated regime is completed, the energy budget reduces in a good approximation to a balance between $E_{p}$ and $E_{d}$ (i.e. $E_{k}\ll E_{p},E_{d}$ ), which is consistent with the assumption that drag and buoyancy forces dominate during the drag-dominated regime. For a given channel Reynolds number, Re, and size of cylinder, $D/H$ (e.g. for $\mathit{Re}=15\,000$ and $D/H=0.032$ ), the temporal variations of the global energy budget terms in (7.2) show that at any given time $E_{k}$ decreases and $E_{p}$ increases with an increase in ${\it\phi}$ (figure 13 a). The number of cylinders per unit surface area increases with ${\it\phi}$ , and as the total drag force and the total dissipation induced by the cylinders are proportional to the number of cylinders per unit area, the current slows down and the rate of dissipation increases.

Figure 13.

The effect of solid volume fraction on the variation of the excess potential energy, $E_{p}$ (dash–dot–dotted line), kinetic energy, $E_{k}$ (solid line), and integral of the total dissipation, $E_{d}$ (dashed line), in the $\mathit{Re}=15\,000$ simulations with ${\it\phi}=0\,\%$ (red line L8), ${\it\phi}=1.25\,\%$ , $D/H=0.032$ (light blue line L2), ${\it\phi}=5\,\%$ , $D/H=0.032$ (blue line L3) and ${\it\phi}=12\,\%$ , $D/H=0.032$ (dark grey line L1). Panel (a) shows the energy variation with time, $t/t_{0}$ , while panel (b) shows the variation of the energy with the front position, $x_{f}/H$ .

Consistent with the results of Ooi et al. (Reference Ooi, Constantinescu and Weber2009) for full-depth currents, the variations of the excess potential energy with time and with the front position are also close to linear during the later stages of the slumping phase in the simulation with ${\it\phi}=0\,\%$ . The decay of the excess potential energy with time is no longer linear in the simulations with ${\it\phi}>1.25\,\%$ . However, when the energy budget terms are plotted as a function of the front position, the variation of $E_{p}$ with $x_{f}/H$ is close to linear and independent of ${\it\phi}$ , provided that the drag induced by the cylinders dominates the basal drag, a regime that is very common for the obstructed flows under consideration. At the end of this section, we use the shallow water model to justify this dependence.

The temporal variation of the dissipation, $E_{d}$ , in the simulations with ${\it\phi}=12\,\%$ and $\mathit{Re}=15\,000$ (figure 14 a) shows that, at any given time, $E_{d}$ is about the same in the simulations with $D/H=0.05$ and 0.07. During the drag-dominated regime, the value of $E_{d}$ is approximately 20 % smaller in the simulation with $D/H=0.032$ compared with the other two simulations. The main reason why at a given time $E_{d}(t)$ is smaller in the simulations with $D/H=0.032$ is because $x_{f}$ is smaller (figure 10 b), so a smaller number of cylinders are situated within the body of the current. However, the variations of $E_{d}$ and $E_{p}$ as a function of $x_{f}/H$ are approximately independent of $D/H$ .

Figure 14.

Temporal variation of the excess potential energy, $E_{p}$ (dash–dot–dotted line), kinetic energy, $E_{k}$ (solid line), and integral of the total dissipation, $E_{d}$ (dashed line): (a) effect of $D/H$ in the $\mathit{Re}=15\,000$ , ${\it\phi}=12\,\%$ simulations with $D/H=0.032$ (dark grey line L1), $D/H=0.05$ (pink line L4) and $D/H=0.07$ (dark blue line L5); (b) effect of Re in the ${\it\phi}=12\,\%$ simulations with $\mathit{Re}=375$ , $D/H=0.05$ (orange line L9), $\mathit{Re}=15\,000$ , $D/H=0.05$ (pink line L4) and $\mathit{Re}=150\,000$ , $D/H=0.05$ (dark green line L10).

In the case of high-Reynolds-number currents advancing in a channel with no obstacles, the front velocity increases monotonically with Re and approaches the theoretical inviscid limit, $U_{f}/u_{b}=0.5$ , at very high Reynolds numbers (see Ooi et al. Reference Ooi, Constantinescu and Weber2009; Tokyay et al. Reference Tokyay, Constantinescu and Meiburg2011). In contrast, in the case of high-Reynolds-number currents propagating in channels containing cylinders, the non-dimensional front velocity at a given time is not a function of $\mathit{Re}$ , provided that the current has undergone a transition to the drag-dominated regime and that the transition to this regime has completed at approximately the same non-dimensional time. For example, this is the case for the simulations ( ${\it\phi}=15\,\%$ , $D/H=0.05$ ) with $\mathit{Re}=150\,000$ and $\mathit{Re}=15\,000$ . The weak influence of Re is because the drag arises mainly from the cylinders and this is close to independent of the Reynolds number, Re, when the cylinder Reynolds number, $\mathit{Re}_{D}$ , is sufficiently high that the drag coefficient, $C_{D}$ , becomes independent of $\mathit{Re}_{D}$ for most of the array.

A justification for the observed nearly linear variation with the front position of $E_{p}/E_{p0}$ (figure 13 b) can be given in terms of the shallow layer formulation

where $L$ is the length of the domain. This implies that $E_{p}/E_{p0}$ must vary linearly with $x_{f}/H$ . Moreover, if $h$ varies linearly with $x$ ( $h/H=(1-x/H)/2$ , which is a very good approximation at least for the low-Reynolds-number cases), then $E_{p}/E_{p0}=1-x_{f}/18H$ . This variation of $E_{p}$ is very close to the numerical result in figure 13(b) for ${\it\phi}>1.25\,\%$ . In the high-Reynolds-number cases, although not exactly linear, $h(x)/H$ remains close to linear (e.g. figure 2 b,c). This is why the variation of $E_{p}/E_{p0}$ with $x_{f}/H$ remains close to linear and varies only weakly with ${\it\phi}$ and Re for the high-Reynolds-number cases with ${\it\phi}>1.25\,\%$ .

8. Summary and conclusions

In this paper, Boussinesq gravity currents initiated from lock-release conditions to propagate though an array of rectangular horizontal cylinders have been investigated using 3D eddy-resolving simulations that compute the flow around individual obstacles and through a two-layer shallow water model. The simulations calculate the evolution of the dependent variables through the current; in particular, this allows the structure of the current, the drag forces and the mixing properties to be evaluated as functions of the non-dimensional parameters describing the array of cylinders. The shallow water description is a simplification of the complete dynamics, but nevertheless it is shown that the bulk dynamics of the currents may be predicted from this framework and that it is often able to reproduce the results of the simulations.

We have discussed the ways in which the size, volume fraction of the cylinders, ${\it\phi}$ , and frontal area of the cylinders per unit volume, $a$ , impact the properties of the flows. Naturally, our detailed analysis only strictly pertains to an array of horizontally aligned square cross-section cylinders that span the flow. However, some of the main trends should be expected to carry over to when the axes of the cylinders are vertical. For example, we anticipate that the coherence of the KH billows will be reduced as ${\it\phi}$ and $aH$ increase, because this depends mainly on the shear stress across the interface, which in turn is dependent upon the magnitude of the drag measured through the combined parameter ${\it\Gamma}_{D}$ . One might expect that the development of the interfacial layer will be similar for different orientations of the cylinders, provided that they remain perpendicular to the streamwise direction.

The use of LES has allowed investigation of scale effects for large-Reynolds-number currents. This is an important contribution, given the scarcity of detailed experimental data on gravity currents propagating in a porous medium at field conditions. That the temporal evolution of the bulk properties of the currents appears to be largely independent of the Reynolds number when $\mathit{Re}>10^{4}$ implies that experimental and numerical studies conducted at moderate scales may accurately represent the behaviour at much larger field scales, even if the details of the vortical structure of the current remain somewhat different.

A main advantage of eddy-resolving simulations that resolve flow past the individual obstacles is that the mean drag coefficient can be directly evaluated once an appropriate definition for the mean velocity within the array has been adopted. The results in the present work showed that if the velocity scale is chosen as the mean velocity within the current from which the wake regions are eliminated, then the predicted values of $C_{D}$ over the quadratic drag-law regime $(1.2<C_{D}<2)$ are fairly comparable with those measured for isolated square cylinders in constant-density uniform flow. Still, for high-Reynolds-number currents these values were found to vary by up to 80 % with the geometrical characteristics of the array ( ${\it\phi},D/H,s/D,aH$ ). Using the approach advocated in the present paper, the drag coefficient can be estimated. This eliminates the need to calibrate simpler numerical models that do not resolve the flow around individual cylinders. Moreover, the predictions of the shallow flow theory model for the temporal evolution of the front position were found to be in good agreement with the simulation results once the drag coefficient was redefined with the mean streamwise velocity within the current. Consistent with shallow water theory, the numerical simulation results confirmed that the combined drag parameter for the porous medium, ${\it\Gamma}_{D}$ , is the key grouping of parameters that determines how fast high-Reynolds-number currents propagate through arrays of obstacles.

An important finding of the present work was that high-Reynolds-number currents propagating through a porous medium are not exactly self-similar; that is, the velocity fields at the front and centre of the flow do not appear to exhibit the same dependences on time. We suggest that a possible explanation for this discrepancy is that the shallow layer model neglects mixing between the dense and less dense fluids. The simulations reveal that the interface is not particularly thick, but that there is significant mixing at the front of the current. Such a process is omitted from the shallow layer framework of § 3 and its inclusion remains an important outstanding challenge.