Lockexchange gravity currents with a low volume of ...

Viewer
Transcript

PUBLICATIONS Journal of Geophysical Research: Oceans RESEARCH ARTICLE 10.1002/2013JC009721 Key Points: Discuss structure of gravity currents in drag-dominated regime Evolution of currents is signiﬁcantly affected by eroded substrate Effect of interfacial waves generated by interaction with eroded substrate

Correspondence to: G. Constantinescu, [email protected] Citation: Tokyay, T., G. Constantinescu, and E. Meiburg (2014), Lock-exchange gravity currents with a low volume of release propagating over an array of obstacles, J. Geophys. Res. Oceans, 119, 2752–2768, doi:10.1002/ 2013JC009721. Received 11 DEC 2013 Accepted 9 APR 2014 Accepted article online 16 APR 2014 Published online 12 MAY 2014

Lock-exchange gravity currents with a low volume of release propagating over an array of obstacles Talia Tokyay1, George Constantinescu1, and Eckart Meiburg2 1 Department of Civil and Environmental Engineering, IIHR-Hydroscience and Engineering, University of Iowa, Iowa City, Iowa, USA, 2Department of Mechanical Engineering, University of California at Santa Barbara, Santa Barbara, California, USA

Abstract The study discusses, based on high-resolution 3-D large eddy simulation results, the evolution of lock-exchange Boussinesq gravity currents with a low volume of release (LVR) propagating over an array of identical two-dimensional (2-D) obstacles (large roughness features in the form of dunes and ribs corresponding to eroded substrates) in a rectangular horizontal channel. The study analyzes the effect of the shape and height of the roughness elements and scale effects between Reynolds numbers at which most laboratory experiments of lock-exchange currents are conducted and Reynolds numbers closer to ﬁeldscale currents in geophysical applications on the temporal variation of the front velocity, mixing, and ﬂow structure within the gravity current. The temporal evolution of the ﬂow instabilities (e.g., Kelvin-Helmholtz billows, forward and backward propagating interfacial waves forming as a result of the interaction of the front of the current with the obstacles) and bed friction velocity distributions are analyzed during the different stages of the evolution of the current. The focus is on the later stages of the propagation of LVR currents, after the transition to the self-similar drag-dominated regime has started. The differences between the evolution and structure of gravity currents with a low and a high volume of release are highlighted.

1. Introduction The study of the propagation of turbulent gravity currents over arrays of obstacles is of practical relevance for many applications in geosciences where, in most cases, gravity currents propagate over rough surfaces [De Cesare et al., 2001; Xu et al., 2004; Meiburg and Kneller, 2010]. This is true for both compositional currents that are driven by differences in salinity and/or temperature between the current and the surrounding ambient ﬂuid and turbidity currents that are driven by suspended sediment. For example, the natural topography over which gravity currents forming in lakes and marine environments propagate contains bed forms (e.g., dunes) and other types of large-scale deformations (e.g., steep-sided macro-roughness observed in the eroded substrate beneath turbidite beds and on the modern see ﬂoor [Eggenhuisen et al., 2010a; Eggenhuisen and McCaffrey, 2012; Bïe et al., 2004; Macdonald et al., 2011]) caused by erosion of the see ﬂoor. Several experiments conducted for constant-ﬂux gravity currents [Sequeiros et al., 2010a, 2010b; Sequeiros, 2012; Eggenhuisen and McCaffrey, 2012] showed that the interaction of the current with the deformed bed or with bottom obstacles results in a signiﬁcant change in their structure with respect to the canonical case of a current propagating over a ﬂat smooth surface corresponding to a plane substrate. One of the most popular ways to control the sediment depositing in reservoirs is the construction of obstacles such as an isolated or a series of submerged dams in the path of the gravity current [Oehy and Schleiss, 2007]. In other applications, man-made-induced gravity currents are used for harbor dredging [Estourgie, 1988]. Most of the information on the structure and evolution of constant-ﬂux and lock-exchange gravity currents comes from scaling relationships, analyses based on simpliﬁed depth-averaged equations [e.g., Garcia, 1993], and laboratory experiments [e.g., Rottman and Simpson, 1983; Stacey and Bowen, 1988; Garcia, 1993; Altinakar et al., 1996; Kneller et al., 1999; Shin et al., 2004; Baas et al., 2005] conducted at Reynolds numbers that are orders of magnitude smaller than those of ﬁeld currents. Most of these experimental investigations were conducted with a plane substrate. Recent laboratory investigations [Sequeiros et al., 2010a; Eggenhuisen and McCaffrey, 2012] considered the case of constant-ﬂux turbidity currents propagating over an eroded bed. Xu et al. [2004] were among the few to provide in situ velocity measurements within turbidity

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2752

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

currents propagating over large bed forms in the Monterey Submarine Canyon and to estimate the average characteristics of the bed forms. Regardless of the shape of the bed, measurements of the instantaneous bed shear stress distributions are nearly impossible to achieve experimentally. Besides laboratory experiments, Reynolds-Averaged Navier-Stokes (RANS) simulations were successfully used to expand the parameter range of laboratory experiments and provide detailed spatial information on the vertical structure of gravity currents [e.g., Huang et al., 2005; Sequeiros, 2012]. Though computationally more expensive, highresolution eddy-resolving numerical simulations [Necker et al., 2002; Ozgokmen et al., 2004; Cantero et al., Figure 1. Sketch of the ﬂow during (a) the slumping phase and (b) the self-similar 2008; Ooi et al., 2009; Tokyay et al., buoyancy-inertia and turbulent drag-dominated phases for the case of a partialdepth lock-exchange ﬂow with a small volume of release (x0/h 5 O(1)). Figure 1c 2012] provide more detailed and shows the shape of the dunes used in the simulations. In some of the simulations, often more accurate information on an array of 2-D ribs or dunes is mounted on the bottom surface to the right of the the turbulent structure of gravity curlock gate. rents. Three-dimensional (3-D) Direct and Large Eddy Simulation based techniques (DNS, LES) have the advantage that they can accurately capture the dynamics of the energetically important eddies and instabilities developing in the ﬂow and thus can be used to better understand the ﬂow physics and to explain trends observed in the mean ﬂow and turbulence statistics (e.g., over regions where the ﬂow within the current is quasisteady [see Tokyay et al., 2012]). In the present study, the evolution of the current and its interaction with the obstacles is studied based on highly resolved 3-D LES performed using a nondissipative Navier-Stokes solver. Such 3-D LES calculations allow simulating currents at Reynolds numbers that are high enough (the Reynolds number deﬁned with the front velocity during the slumping phase and mean height of the head is larger than 10,000) for the ﬂow to be highly turbulent in the head and dissipative wake regions [e.g., Ooi et al., 2007, 2009; GonzalezJuez et al., 2010]. Such currents are referred to as high-Reynolds number currents. This is the case of practical relevance for geophysical applications. The physics of gravity currents propagating over obstacles and large-scale roughness is much more complex than the widely studied case of currents propagating over a ﬂat surface [Rottman et al., 1985; Armi, 1986]. This is true for both constant-ﬂux [e.g., Sequeiros et al., 2010a; Eggenhuisen and McCaffrey, 2012] and lock-exchange gravity currents [e.g., Lane-Serff et al., 1995; Ozgokmen et al., 2004; Ozgokmen and Fisher, 2008; Gonzalez-Juez et al., 2009, 2010; Tokyay et al., 2011a, 2011b, 2012]. Though all lock-exchange currents form as a result of the release of a ﬁnite volume of heavier ﬂuid, in the literature they are often grouped into two categories depending on the initial aspect ratio of the volume of release. In the case of currents with a high volume of release (HVR), the ratio between the initial length and initial height of the release, x0/ h (Figure 1), is much larger than 1. For such cases, after a short acceleration phase, currents propagating over a planar horizontal surface reach the so-called slumping phase, where the front velocity Uf is constant. HVR currents are a convenient way to study currents for which the slumping phase lasts over a long time. Lock-exchange ﬂows for which x0/h is O(1) result in the formation of a so-called current with a low volume of release (LVR). A short time after the formation of the LVR current, the backward propagating current (Figure 1a) hits the endwall of the channel, a reﬂected bore forms and starts advancing toward the front. Once

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2753

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Table 1. Main Parameters of the Simulationsa Case LR-F HR-F LR-R10 LR-R15 HR-R15 LR-D05 LR-D15 HR-D15

Re 47,800 1,000,000 47,800 47,800 1,000,000 47,800 47,800 1,000,000

Obstacle

Rib Rib Rib Dune Dune Dune

xo/(H/2)

h/H

D/(H/2)

k/(H/2)

Position of the Obstacles

0.56 0.56 0.80 0.56 0.56 0.56 0.56 0.56

0.5 0.5 1 0.5 0.5 0.5 0.5 0.5

0.10 0.15 0.15 0.05 0.15 0.15

1 3 3 1 3 3

x/(H/2) 5 3, 4, . . ., 25, 26 x/(H/2) 5 5, 8, 11, 14, 17 x/(H/2) 5 5, 8, 11, 14, 17 x/(H/2) 5 2, 3, ..., 17, 18 x/(H/2) 5 5, 8, 11, 14, 17 x/(H/2) 5 5, 8, 11, 14, 17

a LR denotes Re 5 48,000 simulations while HR denotes Re 5 106 simulations; F, R, and D denote ﬂat bed, ribs, and dunes, respectively; and the number denotes the relative size of the obstacles, D/(H/2).

the bore catches the front, the LVR current transitions to the buoyancy-inertia phase in which Uf decays with time (Figure 1b). In most practical applications (e.g., avalanche-type of currents), x0/h is often very large, but one is interested to analyze the evolution of the current at large times after their formation (e.g., after the end of the slumping phase). This is a main reason why one is interested in studying HVR and LVR currents numerically and in the laboratory. In the case of currents propagating over an array of macro-roughness elements or in a porous medium containing distributed obstacles that induce a large drag force per unit streamwise length compared to the viscous drag at the bed, the evolution of the lock-exchange current is more complex. For example, in the case of high-Reynolds number HVR currents propagating over an array of identical obstacles, the current can transition from the slumping phase to a regime where Uf decays with time. Such a regime in which Uf tb (t is the time measured from the time the lock gate is removed) with b ﬃ 20:28 was observed by Tokyay et al. [2011a] for a HVR current propagating over bottom-mounted square ribs. This regime is similar to the turbulent drag-dominated regime observed for HVR currents propagating into a porous medium for which shallow water ﬂow theory predicts b ﬃ 20:33 [Hatcher et al., 2000]. In this regime, the ﬂow evolution is mainly determined by the balance between the buoyancy and turbulent drag forces on the obstacles within the channel. A drag-dominated regime is also expected to be observed for constant-ﬂux currents subject to a large drag force per unit streamwise length induced by either obstacles present within the channel (e.g., porous or vegetated region) or placed at the channel bottom (e.g., large-scale bed forms). However, the front velocity decay with time is expected to be slower (Uf t21/4 instead of Uf t21/2) for the constantﬂux case [Hatcher et al., 2000]. The obstacles also induce the formation of large-scale instabilities along the interface in the form of backward and forward propagating wave-like disturbances [Armi, 1986; Lane-Serff et al., 1995]. The studies of Tokyay et al. [2011a, 2012] investigated the evolution and structure of HVR currents during the slumping phase and the drag-dominated phase. Ooi et al. [2009] found that important differences were observed between the structure and evolution of LVR and HVR currents propagating over ﬂat smooth surfaces. Some of these differences are expected to be even larger for currents propagating over an array of large-scale obstacles. The present study provides a detailed investigation of the physics of high-Reynolds number compositional LVR gravity currents propagating over an array of 2-D spanwise-oriented identical ﬁxed obstacles (dunes and ribs) on a bed with zero-mean slope. For the reasons discussed before, we focus on the evolution of LVR lock-exchange currents after transition to the buoyancy-inertia phase has started. Though most currents form in regions (e.g., canyons) where the bed slope is high, these currents generally reach a fairly ﬂat region after they move away from the formation region [Garcia, 1993]. In other cases (e.g., lock facilities is an estuarine environment), the ﬁnite-volume current propagates over preexisting bed forms in a region where the average bed slope is generally small. The ﬁndings of this study are strictly valid only for compositional currents. However, the structure of turbidity currents propagating over identical obstacles is expected to be fairly similar to that of the corresponding compositional currents [Sequeiros et al., 2010a]. Tokyay et al. [2012] already discussed some results for LVR currents based on a simulation with large dunes and one with a ﬂat bed (cases LR-D15 and LR-F in Table 1). The discussion was limited to a description of the front velocity variation with time, analysis of the bed friction velocity distributions and sediment entrainment potential of the current propagating over large dunes. In the present paper, we analyze results

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2754

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

of eight simulations (Table 1) with ﬂat beds, dunes, and ribs. This allows us to understand and quantify the effects of obstacle shape, height, and length and Reynolds number-induced effects on the evolution and structure of LVR currents. In particular, in the present study we consider LVR currents whose height near the front remains considerably larger than the obstacle height for large times after their formation, and we show that such currents present important similarities and differences with LVR currents propagating over a ﬂat bed and, respectively, over large obstacles for which the current height is comparable with the obstacle height. Moreover, the present paper presents a detailed analysis of the dynamics and effect of large-scale instabilities on the ﬂow within LVR currents. Simulations of HVR currents were performed by Tokyay et al. [2011a, 2012] for six out of the eight test cases considered in the present study. This allows a clear identiﬁcation of the differences between the physics of LVR and HVR currents propagating over macro-roughness as a function of obstacle shape and of Reynolds number-induced scale effects. The present study does not consider the case of constant-ﬂux currents propagating over an inclined ﬁxed bed, which is very relevant for an important class of geophysical applications. Though the present numerical model can be extended to simulate such cases, such computations are computationally signiﬁcantly more expensive than the case where the bed is parallel to the free surface. Another limitation is that the model does not consider the direct effect of sediment on migration of bed forms. On the other hand, given the extended body of laboratory, theoretical, and numerical work performed for lock-exchange currents propagating over a ﬂat horizontal bed, a detailed study of the effects of an array of identical preexisting bed forms or macro-roughness elements on the structure and evolution of LVR currents is a relevant research topic. Such a study should allow isolating the effect of large-scale roughness on the evolution of LVR currents. The main research questions we try to answer in the present paper are: 1. How does the structure of gravity currents with a ﬁnite volume of release propagating over a rough bed containing macro-roughness elements change during the buoyancy-inertia and drag-dominated regimes? How does the eroding substrate roughness modify the structure of the current compared to the baseline case of a current propagating over a planar substrate (no obstacles)? 2. Supposing the drag force per unit streamwise length is sufﬁciently high such that a high-Reynolds number LVR current transitions to the drag-dominated regime, do we expect that Uf tb with b 5 21/2 independent of the Reynolds number, shape of the obstacles, their spacing, and size relative to the current height? 3. What are the main differences between LVR currents whose height is signiﬁcantly larger than the obstacle height, where the effect of the roughness elements is conﬁned to the lower part of the current, and currents whose height is comparable to that of the obstacles? 4. How do the large-scale instabilities in the ﬂow (e.g., interfacial billows, forward and backward propagating waves generated by the interaction of the current with the obstacles) affect the structure of LVR currents and the bed friction velocity distributions? 5. How does the presence of substrate large-scale roughness affect mixing with respect to the case of a planar substrate? 6. What is the effect of the size and shape of the macro-roughness elements on the bed friction velocity? 7. How do the structure and dynamical properties (e.g., front velocity, propagation of disturbances) of turbulent LVR currents change with the Reynolds number?

2. Description of Numerical Model and Simulations The ﬁnite-volume code [Pierce and Moin, 2001] solves the conservative form of the ﬁltered, dimensionless Navier-Stokes equations (Boussinesq approximation for the density), along with the advection-diffusion equation for the nondimensional concentration, C, on nonuniform Cartesian pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃmeshes. The channel half height, H/2, is used as the length scale. The buoyancy velocity, ub 5 g0 ðH=2Þ, is used as the velocity scale. The reduced gravity is g0 5gðqmax 2qmin Þ=qmax , where qmax and qmin are the initial maximum (lock ﬂuid) and minimum (ambient ﬂuid) densities in the domain and g denotes the gravitational acceleration. The ﬁltered continuity and momentum equations are

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2755

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

@ui 50 @xi @ui @ðui uk Þ @p @ 1 @ui @uk 2Cdi2 52 1 1 1mSGS 1 @xk @xi @xk Re @t @xk @xi

(1) (2)

where p and ui represent the dimensionless pressure and the ﬁltered Cartesian velocity component in the i direction. The coordinates in the three directions are denoted either (x1, x2, x3) in index notation or (x, y, z). The vertical direction is y (i 5 2, in index notation). The variable C in equation (2) is the dimensionless con0 0 0 0 centration C5ðC 0 2C 0min Þ=ðC max 2C min Þ, where C max and C min represent the maximum and minimum initial concentrations in the domain. C is linearly related to the density q. The study considers only Boussinesq currents for which density differences in the ﬂow are small enough such that after linearization the only term containing C in the momentum equation (2) is the gravitational term. C is obtained by solving an advection-diffusion equation @C @ðCuk Þ @ 5 1 @t @xk @xk

1 @C 1aSGS Re Sc @xk

(3)

where Re 5ub ðH=2Þ=m is the Reynolds number and Sc is the molecular Schmidt number deﬁned as the ratio of the molecular viscosity, m, to the molecular diffusivity, j. The subgrid-scale viscosity (mSGS) and subgridscale diffusivity (aSGS) in the subgrid-scale terms in equations (2) and (3) are calculated using a dynamic Smagorinsky model [Pierce and Moin, 2001]. The nondimensional time scale is provided by t0 5 (H/2)/ub. A semi-implicit iterative method that employs a staggered conservative space-time discretization is used to advance the equations in time. A Poisson equation is solved for the pressure. All operators are discretized using central discretizations, except the convective term in the advection-diffusion equation solved for the dimensionless concentration, for which the Quadratic Upwind Interpolation for Convective Kinematics (QUICK) scheme is used. The numerical method is second-order accurate in both space and time. The code was validated for a wide range of constant density and stratiﬁed ﬂows. Chang et al. [2006, 2007] discuss ejection of nonbuoyant and buoyant contaminants from bottom channel cavities. Ooi et al. [2007, 2009] discuss validation for lock-exchange intrusion currents and gravity currents propagating over smooth ﬂat surfaces. Gonzalez-Juez et al. [2009, 2010] analyzed gravity currents impacting isolated circular and square cylinders and successfully validated the numerical predictions of the temporal variation of the drag and lift forces on the obstacle using experimental data. Tokyay et al. [2012] showed that LES predictions of mean velocity and turbulent kinetic energy agreed well with experimental data measured for turbulent gravity currents propagating over a smooth ﬂat surface [Stacey and Bowen, 1988; Kneller et al., 1999]. The Schmidt number is taken to be 600, which corresponds to the saline diffusion in water. Table 1 presents the main parameters of the 3-D simulations discussed in the present study. Simulations conducted at Re 5 47,800 and Re 5 106 are denoted LR and HR, respectively. The presence of a ﬂat bed, dunes, or ribs is indicated by the presence of an F, D, or R in the notation. The nondimensional height of the ribs or dunes, D/(H/2)*100, is indicated by a number. The rib and dune of rank i are denoted as ‘‘rib’’ and ‘‘dune i,’’ respectively. The streamwise positions of the center of the ribs and crest of the dunes along the channel bottom are provided in Table 1. Most of the simulations were conducted with h/(H/2) 5 1 (h is the initial height of the region containing lock ﬂuid) except for case LR-R10 for which h/(H/2) 5 2. No-slip boundary conditions were employed for the velocity at the bottom walls and on the surface of the obstacles. To mimic an open channel environment, the top surface of the channel was treated as a nondeformable slip boundary rather than as a no-slip wall. A convective outﬂow condition [Rodi et al., 2013] which allows the large-scale structures to leave the domain at a ﬁnite speed was used at the right extremity of the computational domain. The surface-normal concentration gradient was set to zero at all wall boundaries. The ﬂow was assumed to be periodic in the spanwise direction, which allows us to study gravity currents of large spanwise extent. The ﬂow ﬁeld was initialized with ﬂuid at rest. The ﬁrst point off the wall surfaces was situated at less than 2 wall units from the bed in the simulations with a ﬂat bed and with ribs. Away from the wall surfaces, the grid spacing in the spanwise and vertical directions was around 0.015(H/2). About 90 points are placed in between the bottom (y/(H/2) 5 0) and y/(H/2) 5 0.2 in the simulations with obstacles of height D 5 0.15(H/2). For Re 5 106, 0.015(H/2) corresponds to about 100 wall units. The time

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2756

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

step was 0.001t0. In most of the simulations, the grid contained 3000 3 64 3 188 nodes in the streamwise (domain length L 5 21.5(H/2)), spanwise (domain width H/2), and vertical (domain height 2(H/2)) directions, respectively. In case LR-R10, the length of the domain was L 5 32(H/2), and the number of grid points in the streamwise direction was close to 4500. Results were checked to be grid independent for the LR-F and LRR15 test cases. The effect of the shape of the obstacles is analyzed by comparing simulations with ribs and dunes of equal height and wavelength (LR-R15 versus LR-D15). The shape of the dunes is taken from the experiment of Mierlo and Ruiter [1988] and is shown in Figure 1c. The k/D ratio of the dunes is about the same as that of the dunes (1–3 m height, 20–70 m wavelength) observed in a ﬁeld study of sandwave migration carried out in the Monterey Canyon by Xu et al. [2004]. In the present study, simulations are conduced with dunes of two different heights, but with constant k/D (520) and of identical shape. Moreover, for the estimated height of the gravity current in this ﬁeld investigation (30 m), the ratio between the obstacle height and the current height was about 0.1, which is the value used in cases LR-D05 and LR-R10 for which D/(h/2) 5 0.1. The case of ribs as macro-roughness elements is also relevant for several reasons. It ﬁrst considers a limiting case in terms of the degree of bluntness of the roughness elements in the array for a given height and wavelength of the elements. Moreover, a bottom-mounted rib of rectangular section was used in a recent laboratory study conducted by Eggenhuisen and McCaffrey [2012] to model a steep-sided scour topography that is often observed at the base of turbidite sandstone beds in the ﬁeld when scouring has removed most of the overlying lithology and left small isolated ridges [Eggenhuisen et al., 2010a, 2010b]. Such positive or negative basal obstructions are also encountered on the sea ﬂoor bathymetry [Wynn et al., 2002]. Finally, arrays of barriers in the form of bafﬂe blocks resembling rectangular ribs are sometimes used to slow down gravity currents in various environments. The effect of the size of the obstacle is mainly studied by comparing results for cases LR-R15 and LR-R10. In the latter case, the region of mixed ﬂuid associated with the gravity current remained signiﬁcantly larger than the height of the ribs during the drag-dominated phase. Additionally, analysis of cases LR-D15 and LRD05 allows comparing an LVR current that transitions rapidly to the drag-dominated regime with a current that does not leave the buoyancy-inertia self-similar phase until the end of the simulation. Scale effects are discussed based on comparison of simulations conducted with Re 5 47,800 (LR-R15) and Re5106 (HR-R15).

3. Temporal Evolution of the Front of the Current Figure 2 shows the temporal evolution of the nondimensional front position (xf-x0)/x0 versus (t/t0) in log-log scale for the low Reynolds number simulations conducted with a ﬂat bed and with ribs of height D 5 0.1(H/2) and D 5 0.15(H/2). As expected, after the end of the initial accelerating phase, the front trajectory in the ﬂat bed cases LR-F (Figure 2) and HR-F (not shown), when plotted in log-log scale, can be well approximated by a straight line of slope equal to one before the reﬂected bore catches the front. The nondimensional front velocity during the slumping phase is Ufs/ub 5 0.53 in case LR-F and 0.58 in case HR-F. In the latter case, the front velocity approaches the value (0.61) predicted by theory [Shin et al., 2004] for partial depth inviscid currents with h/H 5 1/2. The end of the slumping phase corresponds to the location where the straight line with a slope of 1 that approximates well the xf(t) curve plotted in log-log scale over part of the time domain (slumping phase) starts deviating from the curve. Similarly, the start/end of the buoyancy-inertia and drag-dominated self-similar phases correspond to the locations where straight lines with a constant slope of 2/3 and 1/2, respectively, that approximate well xf(t) plotted in log-log scale over part of the time domain start deviating from the xf(t) curve by more than 2%. In some cases of currents propagating over obstacles, the slumping phase is not present and the current transitions directly to the buoyancy-inertia phase (e.g., case LR-R10). During the buoyancy-inertia self-similar phase, simulation results for the ﬂat bed cases show that xf t2/3 (e.g., see Figure 2 for case LR-F), in good agreement with theory [Rottman and Simpson, 1983]. In the simulations with large ribs and dunes, the transition to the buoyancy-inertia self-similar phase ends when (xf-x0)/ x0 5 13–14 in the LR and HR cases. In all these four simulations, the current then transitions to a new regime where the decay of Uf is faster such that xf t1/2. The transition to the new turbulent drag-dominated regime is complete when ðxf -x 0 Þ=x0 ﬃ 20 in the LR cases and when ðxf -x 0 Þ=x0 ﬃ 25 in the HR cases (Table 2). A similar scenario is observed for case LR-R10 where the combined effect of a smaller distance between the obstacles and a smaller height of the obstacles induces a more rapid transition to the

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2757

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

buoyancy-inertia self-similar phase but delays the transition to the drag-dominated regime (ðxf -x 0 Þ=x0 ﬃ 24) compared to the corresponding LR simulations with large obstacles (ðxf -x 0 Þ=x0 ﬃ 20). If the size and degree of bluntness of the obstacles is further reduced (e.g., see results in Table 2 for case LR-D05), the current does not start to transition toward the drag-dominated regime. So for sufﬁciently small-height obstacles and/or obstacles with a large k=D, LVR currents may not transition to a dragdominated regime or this transition will take a very long time. Present results show that high-Reynolds number LVR currents propagating over sufﬁciently large obstacles transition to a dragdominated regime in which xf t1/2 regardless of the shape of the obstacle and Reynolds number. Moreover, our analysis shows that the velocity decay coefﬁcient b 5 20.5 (Uf tb) over the drag-dominated regime predicted by shallow water ﬂow theory and conﬁrmed experimentally for high-Reynolds number LVR currents propagating into a porous medium [Hatcher et al., 2000; Yuksel et al., 2012] is the same for currents propagating over arrays of obstacles, provided that the drag force per (streamwise) unit channel length is sufﬁciently high to induce transition to a turbulent dragdominated regime. However, this coefﬁcient is different than the one observed for HVR currents (c 5 20.28) [Tokyay et al., 2011a, Figure 3].

Figure 2. Time variation of the nondimensional front position, (xf-x0)/x0, as a function of the nondimensional time, t/t0 (log-log scale), in the low Reynolds number simulations with a ﬂat bed (LR-F, red line), with small ribs (LRR10, green line) and large ribs (LR-R15, black line). The dashed lines in the ﬁgure are straight lines with a slope of 1, 2/3, and 1/2 used to identify the slumping, the self-similar buoyancy-inertia, and the drag-dominated phases, respectively.

4. Structure of the Gravity Current 4.1. Currents Propagating Over Large Obstacles In all simulations of LVR gravity currents for which the ratio h/D is small, most of the kinetic energy is concentrated inside the head and dissipative wake after the end of the slumping phase. The presence of large obstacles induces signiﬁcant modulations in the overall decay of the velocity with the distance from the front. If the obstacles are sufﬁciently large, similar to the case of HVR currents [see Tokyay et al., 2011a], a jet-like ﬂow forms past the top of the obstacles situated close to the front because of the strong acceleration of the heavier ﬂuid convected over the top of the obstacle as it approaches the lower parts of the bed. The additional drag induced by the obstacles reduces the velocity over the whole length of the current and the mixing induced by the interaction of the front with the obstacles decreases the mean value of the concentration/density in the head region compared to the case of an LVR current propagating over a ﬂat surface. The energy associated with the large-scale eddies generated by the interaction of the current with the

Table 2. Comparison of the Nondimensional Times and Front Positions Corresponding to the End of the Slumping Phase, the Start of the Buoyancy-Inertia Self-Similar Phase, and the Start of the Drag-Dominated Self-Similar Phase End of Slumping Phase

LR-F HR-F LR-R10 LR-R15 HR-R15 LR-D05 LR-D15 HR-D15

TOKYAY ET AL.

Start of Buoyancy-Inertia Phase

xf 2xo xo

xf H=2

t to

xf 2xo xo

xf H=2

11 9

9.5 9.6

5.9 5.9

9 7.5 10 9 7.5

7.8 8 8.5 7.7 8

4.9 5 5.3 4.9 5

25 25 10 16.5 15.5 22 17.3 16.6

18.3 22 8.2 12.4 13.7 16.3 13.6 14.8

10.8 12.9 7.4 7.5 8.2 9.7 8.2 8.8

t to

C 2014. American Geophysical Union. All Rights Reserved. V

Start of Drag-Dominated Phase t to

xf 2xo xo

xf H=2

51 30 42

24 18.5 24.9

20 11 14.5

32 40

21 25.8

12.3 15

2758

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Figure 3. Structure of the current in the case LR-R15 during the drag-dominated phase when xf =ðH=2Þ ﬃ 16. (top) The distributions of the concentration in a x-y section. (bottom) The streamwise distribution of the Froude number.

obstacles and their capacity to mix vertically ﬂuid within the current increases with the relative size of the obstacles. In this regard, the role of the large-scale eddies induced by the interaction of a current with a bottom obstacle is similar for the LVR currents considered in this study and for constant-ﬂux currents interacting with an isolated rectangular 2-D rib placed at the channel bed [Eggenhuisen and McCaffrey, 2012]. 1=2 The streamwise distributions of the Froude number, Fr5ð u =ub Þ=ðh=ðH=2ÞÞ , and current height, h=ðH=2Þ, are given in Figures 3–6. The current height and the mean streamwise velocity at a certain streamwise locaÐH Ð Ð =ub 5 0H ðu=ub ÞCdy= 0H Cdy, respectively. tion are estimated as h=ðH=2Þ51=ðH=2Þ 0 Cdy and u

In the simulations with large obstacles (D50.15(H/2)), the Froude number in the head region (e.g., see Figure 3 for case LR-R15) is very close to unity during the buoyancy-inertia and drag-dominated self-similar phases. The ﬂow strongly accelerates within the jet-like ﬂow region forming past the top/crest of the ﬁrst 1–2 obstacles behind the front. This is similar to what was observed for HVR currents by Tokyay et al. [2011a, 2012]. The intensity of the jet-like ﬂow decreases sharply for the obstacles situated behind the head once the current transitions to the drag-dominated regime. The trapping of higher-density ﬂuid in between successive obstacles explains why the average height of the current between two obstacles decreases with the rank of the obstacle in Figures 4b and 4c. By contrast, in the ﬂat bed cases (Figure 4a) the height of the current increases slightly with x in between the rear endwall and the back of the dissipative wake. The mixing induced by the interaction of the front with the obstacles causes a signiﬁcant decrease in the average heights of the head and dissipative wake in the simulations with obstacles compared to the corresponding ﬂat bed cases (e.g., by about 40% between cases LR-R15 and LR-F and by about 60% between cases LR-R15 and LR-F when xf =ðH=2Þ ﬃ 18, as shown in Figure 4).

Figure 4. Streamwise variation of the nondimensional height of the gravity current, h=ðH=2Þ, when xf =ðH=2Þ ﬃ 23:5. (a) Cases LR-F (solid line) and HR-F (dashed line); (b) cases LR-R15 (solid line) and HR-R15 (dashed line); (c) cases LR-D15 (solid line) and HR-D15 (dashed line).

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2759

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Figure 5. Structure of the current in the case LR-R10. (a) Distribution of the concentration in a x-y plane during the buoyancy-inertia selfsimilar phase when xf =ðH=2Þ ﬃ 11; (b–d) distribution of the concentration and streamwise velocity in a x-y plane and of the Froude number during the drag-dominated phase when xf =ðH=2Þ ﬃ 22:5.

The main effect of augmenting the Reynolds number from 48,000 to 106 is an increase of the height of the current in the head and dissipative wake regions (Figure 4a). Relatively speaking, the increase is larger in the case of a ﬂat bed and for the blunt obstacles (ribs). This happens because mixing close to the front is signiﬁcantly less in the higher-Reynolds number case. In the case of ribs, the amount of lower density ﬂuid trapped beneath the higher-density ﬂuid originating in the splash forming as the front interacts with a new rib decreases signiﬁcantly with the increase of the Reynolds number. The increased energy of the highly 3D eddies within the head and dissipative wake regions in the higher-Reynolds number simulations reduces the size and coherence of the KH billows shed behind the head compared to the corresponding lower Reynolds number simulations with a ﬂat bed and with obstacles. In the simulations with large obstacles, the effect of the Reynolds number on the height variation in between the crest/top of two consecutive obstacles situated behind the front is negligible (Figures 4b and 4c). This happens because the concentration distribution in between successive obstacles situated at large distances behind the front is mainly a function of the amount of higher-density ﬂuid trapped in between the top/crest of the obstacles which, for sufﬁciently large obstacles, is primarily a function of the shape and height of the obstacles.

Figure 6. Streamwise variation of the nondimensional height of the gravity current, h=ðH=2Þ, in the case LR-R10. (a) xf =ðH=2Þ ﬃ 10:5; (b) xf =ðH=2Þ ﬃ 23:5.

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2760

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

4.2. Currents Propagating Over Small Obstacles Analysis of the structure of the current in case LR-R10 is very signiﬁcant for cases in which the height of the region containing ﬂuid of density larger than that of the ambient ﬂuid close to the front remains larger than the obstacle height even after the end of the transition to the drag-dominated regime. The value of h/D (510) in case LR-R10 is much larger than the one in the other simulations with obstacles (e.g., LR-R15) for which h/D 5 3.33. No jet-like ﬂow was detected on the lee side of the ribs in case LR-R10. The height of the region containing mixed ﬂuid remains much larger than the size of the ribs during most of the buoyancyinertia phase. The shape of the region containing the head and dissipative wake resembles the one observed for currents propagating over a smooth ﬂat surface (Figure 5a). Most of the differences with a LVR current propagating over a smooth surface are conﬁned to the near bed region. Gradually, due to successive mixing induced by the ribs and the reduction in the height of the front region, the ﬂuid inside the head gets more and more diluted. For the conditions considered in case LR-R10, the concentration/density levels in the region of high streamwise velocities associated with the head become much lower than those inside the tail, after the start of the transition to the drag-dominated regime (Figure 5b). The presence of closely spaced ribs and the absence of a jet-like ﬂow past the ribs cause the region of high streamwise velocity within the current to be situated mostly above the level corresponding to the top of the ribs, with little penetration of high-velocity ﬂuid near the bed, in between the ribs (e.g., see Figure 5c). The distribution of the current height before the start of the transition to the drag-dominated regime (Figure 6a) resembles more the one observed for LVR currents propagating over a smooth surface during the buoyancy-inertia self-similar phase (Figure 4a) for which h=ðH=2Þ peaks over the head and dissipative wake. During the drag-dominated regime a close to linear decrease of the Froude number behind the front of the current is observed. For example, in Figure 5d the linear decay regime ends at x/(H/2) 5 12 where Fr ﬃ 0.

5. Dynamics of Flow Disturbances and Interfacial Billows The x-t plots of the concentration C on a line situated at a small distance (y/(H/2) 5 0.09) from the channel bottom in Figure 7 visualize the front trajectory and the main forward and backward propagating disturbances in the ﬂat bed case (LR-F) and in the case with small dunes (LR-D05). The front trajectory corresponds to the boundary between the white region containing only ambient ﬂuid and the darker region. The slope is approximately constant up to x/(H/2) 5 5.3–6 when the reﬂected disturbance overtakes the front (Table 2). As Uf 5 dxf/dt starts decaying, the front trajectory curves toward the time axis. The dark streaks present close to the front in Figure 7 visualize the trajectories of the large interfacial billows being shed behind the head. These billows have initially a velocity close to that of the front but then gradually reduce their velocity. Before they dissipate, the angle of the streaks with the time axis becomes quite small, indicating that these billows move downstream with a very small velocity at large times after their formation. Figure 7 also shows the presence of three streaks of darker color denoted BD1, FD1, and FD2. They correspond to the trajectories of several disturbances that propagate within the gravity current (see also Figure 13a). BD1 and FD2 originate at a common point in Figure 7 situated at (t ﬃ 4t 0 , x ﬃ H=2). A dark streak corresponding to a large billow separating from the dissipative wake passes through that point in the C(x,t) plot. Animations show that the core of that billow starts moving toward the bed as it is advected downstream with a very small velocity. As it interacts with the channel bottom, some of the higher-density ﬂuid contained within the billow starts moving downstream creating the forward-moving disturbance front FD2. The other part of the higher-density ﬂuid starts moving toward the rear endwall in the form of a backwardmoving disturbance front BD1. As BD1 approaches the rear endwall ðt=t0 ﬃ 10Þ, it reﬂects as a forwardmoving disturbance, FD1. These propagating waves give rise to local perturbations of the current depth and concentration. The disturbances transfer energy and momentum along the interface. The trajectories of BD1, FD1, and FD2 are close to straight lines. For case LR-F, the corresponding velocities of the BD1, FD1, and FD2 are 0.14ub, 0.145ub, and 0.25ub, respectively. Though the front velocity is higher in case HR-F, the disturbances propagate at slightly lower speeds (0.13ub, 0.14ub, and 0.21ub). The presence of small obstacles with a reduced degree of bluntness in case LR-D05 has only a negligible effect on the propagation of the ﬂow disturbances compared to the ﬂat bed case, LR-F (Figure 7). The forward-propagating disturbances transfer energy toward the current front. The pattern of higher-concentration streaks is much more complex in the simulations with obstacles (e.g., see Figure 8 for case LR-R15). Similar to case LR-F, a large billow containing higher-density (lock) ﬂuid

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2761

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Figure 7. Evolution of the concentration with time along the x axis. (a) Case LR-F; (b) case LR-D05. The distance from the bottom wall is y/ (H/2) 5 0.09.

separates from the dissipative wake region. As the billow reaches the channel bottom around x 5 H/2, it gives rise to a stronger forward-propagating disturbance, S1, and to a weaker backward-propagating disturbance, BD1. BD1 reaches the rear endwall when t 12t0 and then reﬂects as a forward-propagating disturbance, FD12 (Figure 9). Another forward-propagating disturbance, FD11, forms when BD1 encounters a patch of higher-density ﬂuid that separated earlier from the back of the current and that moves downstream. The strengths of both FD11 and FD12 are much smaller than those of S1 and of the backwardpropagating disturbance S2 that forms as the front starts interacting with rib1 at t 10t0. These disturbances are visualized in Figure 9. Each time the front reaches a rib, it gives rise to a backward-propagating disturbance (S3 for rib2, S4 for rib3 and S5 for rib4). The velocity of the disturbance S1 as it propagates toward rib1 is close to 0.25ub, which is about half that of the current front during the slumping phase. The velocity of S1 after reﬂection at rib1 is about 0.24ub. The velocity of S2 as it propagates from rib1 toward the endwall is slightly higher (0.27ub). The subsequent reﬂections of S1 and S2 at the rear endwall are accompanied by a noticeable loss in the disturbance velocity (e.g., from 0.27ub to 0.19ub for S2). Analysis of the concentration and streamwise velocity ﬁelds allows understanding how two disturbances propagating in opposite directions pass each other. For example, in Figure 10a, S1 moves downstream while S3 moves upstream. As the disturbances approach each other (Figure 10b), S3, that contains lower density ﬂuid compared to S1, moves over S1. After S1 has passed over S3, its center approaches again the channel bottom, as shown by the position of the region of high positive streamwise velocity associated with S1 in Figure 10c. After they intersect, the velocity of S1 decreases quite fast, while that of S3 is left unchanged (Figure 8). Figure 8. Evolution of the concentration with time along the x axis in the case LR-R15. The distance from the bottom wall is y/(H/2) 5 0.09. The dashed lines visualize the main forward and backward propagating ﬂow disturbances. S1 and BD1 form due to the collapse of a large interfacial billow on the channel bottom. S2–S5 form due to the partial reﬂection of the gravity current ﬂuid as the front passes Rib1–Rib4, respectively.

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

Figures 10 and 11 give more details on the dynamics of the disturbances and their effect on the structure of the gravity current, as they interact with rib2. The dynamics of the

2762

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Figure 9. Visualization of the structure of the current in between the rear endwall and rib1 at t 5 21t0 for case LR-R15 using (left) nondimensional concentration contours and (right) streamwise velocity contours. Also shown are the position and direction of propagation of the main ﬂow disturbances.

disturbances is qualitatively similar for the higher-ranked ribs. As the backward-propagating disturbance reaches the rib situated upstream of it, the disturbance reﬂects and starts moving forward with a smaller velocity (see length of time intervals needed for each disturbance to travel between the same two ribs in Figure 8). For example this is the case for S3. Figures 10b and 10c show S3 before and after reﬂection at rib1. As S3 reaches rib1, a small amount of higher-concentration ﬂuid is pushed over rib1 and continues to propagate toward the endwall. The other part gets reﬂected and starts propagating toward rib2. Another interesting situation occurs around t 5 58t0 when the forward-propagating disturbance S3 and the backward-propagating disturbance S4 reach rib2 at about the same time (Figure 11). Both disturbances have enough momentum to carry higher-density ﬂuid from the vicinity of the rib past its top. The net amount of higher-density ﬂuid convected over the top of the rib is very small. The trajectories of the main disturbances S1–S5 are similar in cases LR-R15 and HR-R15.

6. Mixing The presence of bottom obstacles inﬂuences mixing between the bottom current and the ambient ﬂuid. Information on mixing is important for geophysical applications. For example, LVR currents which mix faster with the ambient ﬂuid will loose more rapidly their capacity to entrain sediment from the loose bed in the later stages of their propagation. A simple way to quantify mixing [e.g., Necker et al., 2005] is to estimate the variation with time of the subvolume containing ﬂuid with density/concentration higher than that of the lighter ambient ﬂuid. After nondimensionalization with the initial volume of lock ﬂuid, V0, the normalized subvolume containing ﬂuid with a density that is higher than that of the ambient ﬂuid can be calculated as Ð V’ðxf =HÞ5 V10 X adV, where a 5 1 if 0.002 < C < 1 and a 5 0 otherwise. Initially, the volume of the LVR current is V0, so V0(0) 5 1. In the case of large bottom obstacles, mixing is larger for a LVR current propagating over ribs compared to the case when the same current propagates over dunes. For example, in the case of obstacles with a height D 5 0.15(H/2), results in Figure 12 show that V0 is by about 5% larger in the simulation with ribs after the transition to the drag-dominated regime has started. In the case when obstacles with a high degree of

Figure 10. Visualization of the structure of the current around rib1 and rib2 for case LR-R15 using (left) nondimensional concentration contours and (right) streamwise velocity contours. Also shown are the position and direction of propagation of the main ﬂow disturbances. (a) t 5 26t0; (b) t 5 30t0; (c) t 5 39t0.

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2763

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Figure 11. Visualization of the structure of the current around rib2 for case LR-R15 using (left) nondimensional concentration contours and (right) streamwise velocity contours. Also shown are the position and direction of propagation of the main ﬂow disturbances. (a) t 5 52t0; (b) t 5 560; (c) t 5 62t0.

bluntness (ribs) are present at the channel bed, the splash of mixed ﬂuid generated by the interaction of the front with the upstream face of the obstacle and the trapping of ambient ﬂuid behind the rib after the plunging head reattaches to the channel bed explain the larger values of V0 in case LR-R15 compared to case LR-F. In the case when large obstacles with a low degree of bluntness (dunes) are present at the channel bed, one of the mechanisms which promotes mixing in the simulation with ribs is absent, as the interaction of the front with the upstream face of the dune does not generate a splash. However, a larger amount of ambient ﬂuid is trapped in between the nose of the current and the bed surface when the front advances over the downslope part of each dune compared to the case when the current propagates over a ﬂat surface. Moreover, the propagation of the disturbances generated by the interaction of the front with the obstacles can also enhance mixing around the interface of the current. In the case of an LVR current propagating over small obstacles, the effect of the bottom roughness on the total volume of mixed ﬂuid is relatively small. This is especially the case for obstacles with a low degree of bluntness. For example, the variation of V0 for case LR-D05 (not shown) was found to be very close to the one observed for the ﬂat bed case (LR-F in Figure 12).

7. Bed Friction Velocity

Figure 12. Mixing of the lock ﬂuid and ambient ﬂuid in the LR-F, LR-D15, and LR-R15 simulations. The mixing is quantiﬁed by the normalized volume of mixed ﬂuid V0 (0.002 < C < 1) as a function of the front position, xf/H.

TOKYAY ET AL.

Similar to HVR currents [Tokyay et al., 2012], the main region of high us is situated behind the front of LVR currents propagating over a ﬂat surface (Figure 13a). Streaks of low and high us are observed for some distance behind the front in the region where the ﬂow is highly turbulent. Moreover, as shown by the 1-D proﬁles of the spanwise-averaged bed friction velocity magni~ s , for case LR-F in Figure 15a, the decrease tude, u ~ s over the head and dissipative wake region is of u close to linear after the start of the buoyancyinertia self-similar phase. The deviations from the linear decay are induced by the interfacial billows shed from the dissipative wake and by the passage of ﬂow disturbances. The other interesting feature of the distribution of us at large distances behind the head in Figure 13a is the formation of a spanwise band of relatively high us between x/ (H/2) 5 1 and x/(H/2) 5 1.8 due to the passage of

C 2014. American Geophysical Union. All Rights Reserved. V

2764

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Figure 13. Spatial distribution of the bed friction velocity magnitude, us/ub, during the buoyancy-inertia self-similar phase when xf =ðH=2Þ ﬃ 9. (a) Case LR-F; (b) case LR-R15. Also shown are the concentration contours in an x-y section.

the forward-propagating disturbance FD1 (Figure 7a). The local increase of the near-bed streamwise velocity in this region is sufﬁciently large to generate streaks of low and high us, similar to those observed behind the front, in a region where the ﬂow started to relaminarize. All these features of the us distributions are also observed for case HR-F (not shown). The effect of the obstacles for LVR currents with small h/D (case LR-R15) and large h/D (case LR-R10) is discussed next. If the spacing between the ribs is relatively large (case LR-R15), a region of high us containing streaks is generally observed behind the front except during and some time after the front passes a new rib (e.g., streaks are not present behind the front in Figure 13b). Large values of us, which can sometimes be comparable to those behind the front, are also observed in the region where the jet-like ﬂow past the ﬁrst 1–2 ribs behind the front reaches the bed surface. Finally, similar to the ﬂat bed cases, the passage of disturbances can strongly increase the near-bed velocity. For example, the passage of the backward-propagating disturbance S2 (Figure 8) explains the presence in Figure 13b of a region of high us containing streaks around x/(H/2) 5 1.8. The streamwise velocity is negative in the immediate vicinity of the bed, even though the overall movement of the ﬂuid contained within the tail in the surrounding region is in the positive streamwise direction. The near-bed ﬂow is strongly turbulent over the head and dissipative wake in case LR-R10, as shown by the high variability of us in the spanwise direction within these regions (Figure 14). The level of turbulence in the near-bed region is higher in case LR-R10 compared to case LR-R15, mainly because of the small, but very energetic, turbulent eddies forming in the shear layers past the ribs. The passage of ﬂow disturbances results in a local ampliﬁcation of us. However, the smaller size of the obstacles induces the formation of much weaker backward-propagating disturbances compared to case LR-R15. If one discards the local varia~ s =ub around each rib and the temporary ampliﬁcation associated with the passage of a ﬂow distion of u turbance, the results in Figure 15a show that, after the end of the transition to the buoyancy-inertia selfsimilar phase, the decay of the bed friction velocity in case LR-R10 is close to linear over the head, dissipative wake and the downstream part of the tail region. So in this regard, gravity currents propagating over obstacles with large h/D are similar to gravity currents propagating over a ﬂat bed [Tokyay et al., 2012]. ~ s =ub is comparable in the simulations The effect of varying the size of the obstacle on the distributions of u with dunes [Tokyay et al., 2012, Figures 7b and 9] and with ribs (Figure 15). For sufﬁciently small dunes (case ~ s =ub with the distance from the front is close to linear beneath the head and dissipaLR-D05), the decay of u tive wake regions. This is similar to what was observed for case LR-R10. ~ s , in For a given front position, the spanwise-averaged distributions of the bed friction velocity magnitude, u the low and high-Reynolds number simulations can be considered to be close after rescaling with a constant factor, c (Figure 16). This result is very relevant for applications where one needs to estimate the sediment entrainment capacity of gravity currents at ﬁeld conditions. Data obtained from numerical simulations and laboratory experiments conducted for turbulent currents at lower Reynolds numbers can be used to

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2765

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Figure 14. Spatial distribution of the bed friction velocity magnitude, us/ub, in the case LR-R10 during the ﬁnal stages of the transition to the drag-dominated phase when xf/(H/2) 5 16.5. Also shown are the concentration contours in an x-y section.

estimate this quantity. The rescaling is needed because the velocity gradients at the bed increase with Re, while us scales with the square root of the molecular viscosity, or Re21/2 when the molecular viscosity is expressed in nondimensional form. Similar to the case of HVR currents [Tokyay et al., 2012], the same value ~ s ðxÞ for a given front position in the simulations with a ﬂat bed, with ribs of c (0.43) can be used to rescale u and with dunes. Tokyay et al. [2012] proposed an approximate way of estimating the value of c for HVR currents based on the value of the bed friction velocity estimated using the Moody diagram for fully developed ﬂow over a smooth bed. A physical Reynolds number with the front velocity and the height of the head is ~ s decays linearly with the distance used to estimate the friction factor. Similar to case LR-F, a region where u from the front is also present in case HR-F. For the ﬂat bed cases, the length of this region increases with the Reynolds number (Figure 16a).

8. Summary and Conclusions High-resolution 3-D LES simulations allowed clarifying some important aspects related to the propagation of Boussinesq turbulent gravity currents with a low volume of release (LVR) over smooth ﬂat surfaces (planar substrate) and over surfaces containing an array of identical obstacles (rough substrate containing macroroughness) in the form of 2-D ribs and dunes. In particular, the present simulations demonstrated the important role played by the backward-propagating disturbances forming at the crest/top of the obstacles and by their reﬂections in transferring momentum and energy between the upstream and the downstream part of the current and in locally increasing the bed friction velocity. Results showed that high-Reynolds number LVR gravity currents propagating over a series of obstacles may transition to a drag-dominated phase in which the front velocity decays proportional to ta, with a 5 21/2. Provided that the additional form drag induced by the obstacles is much larger than the friction drag, the value of a was shown to be independent of the shape and size of the obstacles and of the Reynolds number. For currents propagating over macro-roughness elements, there is no clear way to deﬁne the equivalent volume fraction accounting for the presence of the large-scale roughness elements in the bottom layer. This makes difﬁcult to propose analytical models similar to the ones available for LVR currents

Figure 15. Streamwise distribution of the spanwise-averaged bed friction velocity magnitude in the cases LR-F (red line), LR-R15 (green line), and LR-R10 (black line). (a) xf/(H/2) 5 13.5–14.5; (b) xf/(H/2) 5 23.5. The black and red dashed lines show the region where the overall decay of the bed friction velocity in the cases LR-R10 and LR-F is close to linear. The solid vertical lines show the position of the ribs in the case LR-R15. The vertical arrows show the position of the front.

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2766

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

Figure 16. Streamwise distribution of the spanwise-averaged bed friction velocity magnitude. (a) Cases LR-F (solid line) and HR-F (dashed line); (b) cases LR-R15 (solid line) and HR-R15. The values of us in the low Reynolds number simulations (cases LR-F and LR-R15) are multiplied by a factor c 5 0.43. The vertical arrow shows the position of the front.

advancing into a porous medium [Hatcher et al., 2000] that can be used to estimate the times at which the current transitions to the buoyancy-inertia self-similar phase and to the drag-dominated phase, as well as the shape of the current during these regimes. High-resolution simulations can provide such information. Similar to the case of HVR currents, the regions of high bed friction velocity magnitude for LVR currents with a low h/D corresponded to the ones where the jet-like ﬂow forming downstream of the top/crest of the obstacle plunged toward the lower regions of the bed and to a limited distance behind the front. As opposed to HVR currents, large values of us were induced downstream of only the ﬁrst couple of obstacles behind the front. The evolution and structure of LVR currents with a high h/D ratio were found to share similarities with both LVR currents propagating over a ﬂat smooth surface and LVR currents with a low h/D. Similar to LVD currents propagating over a ﬂat surface during the buoyancy-inertia self-similar phase, the bed-friction velocity decay was close to linear in between the front region and the location where the streamwise velocity in the tail becomes negligible. Moreover, for LVR currents with a high h/D, the streamwise decay of the local Froude number was also close to linear over the same region. Regardless of the value of h/D, the collapse of some of the large Kelvin-Helmholtz billows shed at the back of the LVR current on the channel bottom during the initial stages of the propagation of the current was found to induce the formation of strong forward and backward instabilities. The effect of increasing the Reynolds number in the simulations of LVR currents propagating over relatively large roughness elements was to increase the nondimensional distance traveled by the front before the transition to the drag-dominated regime starts. The speeds of the disturbances induced by the interaction of the front with the roughness elements were fairly insensitive to the Reynolds number in the simulations with a rough substrate containing relatively large macro-roughness elements. For a given position of the front, the mean concentration (excess density) level in the head of LVR currents decreased with the increase in the Reynolds number. However, Reynolds number effects on the variation of the height of the current behind the head region were not signiﬁcant. With the exception of inducing a decay by a close to constant ratio of the overall levels of us with the Reynolds number, no important scale effects were observed between the low (Re 5 48,000) and the high (Re 5 106) Reynolds number simulations. This was true for both the case of a ﬂat substrate and that of a rough substrate. Acknowledgments We gratefully acknowledge the National Center for High Performance Computing (NCHC) in Taiwan, in particular W. F. Tsai, and the Transportation Research and Analysis Computing Center (TRACC) at the Argonne National Laboratory, in particular S. Lottes, for providing substantial computer time. Eckart Meiburg acknowledges ﬁnancial support through NSF grants CBET0854338, CBET-1067847, and OCE1061300.

TOKYAY ET AL.

References Altinakar, M. S., W. H. Graf, and E. J. Hopﬁnger (1996), Flow structure in turbidity currents, J. Hydraul. Res., 34, 713–718. Armi, L. (1986), The hydraulics of two ﬂowing layers with different densities, J. Fluid Mech., 163, 27–58. Baas, J. H., W. D. McCaffrey, P. D. W. Haughton, and C. Choux (2005), Coupling between suspended sediment distribution and turbulence structure in a laboratory turbidity current, J. Geophys. Res., 110, C11015, doi:10.1029/2004JC002668. Bïe, R., T. Bugge, L. Rise, G. Eidnes, A. Eide, and E. Mauring (2004), Erosional channel incision and the origin of large sediment waves in Trondheimsfjorden, Central Norway, Geo Mar. Lett., 24, 225–240. Cantero, M., S. Balachandar, M. Garcia, and D. Bock (2008), Turbulent structures in planar gravity currents and their inﬂuence on ﬂow dynamics, J. Geophys. Res., 113, C08018, doi:10.1029/2007JC004645. Chang, K., G. Constantinescu, and S. O. Park (2007), The purging of a neutrally buoyant or a dense miscible contaminant from a rectangular cavity. Part II: The case of an incoming fully turbulent overﬂow, J. Hydraul. Eng., 133(4), 373–385. Chang, K. S., G. Constantinescu, and S.-O. Park (2006), Analysis of the ﬂow and mass transfer processes for the incompressible ﬂow past an open cavity with a laminar and a fully turbulent incoming boundary layer, J. Fluid Mech., 561, 113–145.

C 2014. American Geophysical Union. All Rights Reserved. V

2767

Journal of Geophysical Research: Oceans

10.1002/2013JC009721

De Cesare, G., A. Schleiss, and F. Hermann (2001), Impact of turbidity currents on reservoir sedimentation, J. Hydraul. Eng., 127(1), 6–16. Eggenhuisen, J. T., and W. D. McCaffrey (2012), The vertical turbulence structure of experimental turbidity currents encountering basal obstructions: Implications for vertical suspended sediment distribution in non-equilibrium currents, Sedimentology, 59, 1101–1120. Eggenhuisen, J. T., W. D. McCaffrey, P. D. W. Haughton, and R. W. H. Butler (2010a), Small-scale spatial variability in turbidity current ﬂow controlled by roughness resulting from substrate erosion: Field evidence for a feedback mechanism, J. Sediment. Res., 80, 129–136. Eggenhuisen, J. T., W. D. McCaffrey, R. W. H. Butler, and P. D. W. Haughton (2010b), Shallow erosion beneath turbidity currents and its impact on the architectural development of turbidite sheet systems, Sedimentology, 58, 936–959. Estourgie, A. L. P. (1988), Theory and practice of water injection dredging, Terra Aqua, 38, 21–28. Garcia, M. H. (1993), Hydraulic jumps in sediment-driven bottom currents, J. Hydraul. Eng., 119, 1094–1117. Gonzalez-Juez, E., E. Meiburg, and G. Constantinescu (2009), Gravity currents impinging on bottom mounted square cylinders: Flow ﬁelds and associated forces, J. Fluid Mech., 631, 65–102. Gonzalez-Juez, E., E. Meiburg, T. Tokyay, and G. Constantinescu (2010), Gravity current ﬂow past a circular cylinder: Forces and wall shear stresses and implications for scour, J. Fluid Mech., 649, 69–102. Hacker, J., P. F. Linden, and S. B. Dalziel (1996), Mixing in lock-release gravity currents, Dyn. Atmos. Oceans, 24, 183–195. Hatcher, L., A. J. Hogg, and A. W. Woods (2000), The effects of drag on turbulent gravity currents, J. Fluid Mech., 416, 297–314. Huang, H., J. Imran, and C. Pirmez (2005), Numerical model of turbidity currents with a deforming bottom boundary, J. Hydraul. Eng., 131, 283–293. Kneller, B., S. J. Bennett, and W. D. McCaffrey (1999), Velocity structure, turbulence and ﬂuid stresses in experimental gravity currents, J. Geophys. Res., 104, 5381–5391. Lane-Serff, G. F., L. M. Beal, and T. D. Hadﬁeld (1995), Gravity current ﬂow over obstacles, J. Fluid Mech., 292, 39–53. Macdonald, H. A., R. B. Wynn, V. A. I. Huvenne, J. Peakall, D. G. Masson, P. P. E. Weaver, and S. D. McPhail (2011), New insights into the morphology, ﬁll, and remarkable longevity (>0.2 m.y.) of modern deep-water erosional scours along the northeast Atlantic margin, Geosphere, 7, 845–867. Meiburg, E., and B. Kneller (2010), Turbidity currents and their deposits, Annu. Rev. Fluid Mech., 42, 135–156. Mierlo, M. C., and J. C. de Ruiter (1988), Turbulence measurements over artiﬁcial dunes, Rep. Q789, Delft Hydraul. Lab., Delft, Netherlands. Necker, F., C. H€artel, L. Kleiser, and E. Meiburg (2002), High-resolution simulations of particle-drive gravity currents, Int. J. Multiphase Flow, 28, 279–300. Oehy, C., and A. Schleiss (2007), Control of turbidity currents in reservoirs by solid and permeable obstacles, J. Hydraul. Eng., 133, 637–648. Ooi, S. K., S. G. Constantinescu, and L. Weber (2007), A numerical study of intrusive compositional gravity currents, Phys. Fluids, 19, 076602, doi:10.1063/1.2750672. Ooi, S. K., S. G. Constantinescu, and L. Weber (2009), Numerical simulations of lock exchange compositional gravity currents, J. Fluid Mech., 635, 361–388. Ozgokmen, T. M., and P. F. Fisher (2008), On the role of bottom roughness in overﬂows, Ocean Modell., 20, 336–361. Ozgokmen, T. M., P. F. Fischer, J. Duan, and T. Iliescu (2004), Entrainment in bottom gravity currents over complex topography from threedimensional nonhydrostatic simulations, Geophys. Res. Lett., 31, L13212, doi:10.1029/2004GL0200186. Pierce, C. D., and P. Moin (2001), Progress-variable approach for large-eddy simulation of turbulent combustion, Mech. Eng. Dep. Rep. TF-80, Stanford Univ., Palo Alto, Calif. Rodi, W., G. Constantinescu, and T. Stoesser (2013), Large Eddy simulation in hydraulics, IAHR Monograph, CRC Press, Taylor & Francis Group (ISBN-10: 1138000247). Rottman, J. W., and J. E. Simpson (1983), Gravity currents produced by instantaneous releases of a heavy ﬂuid in a rectangular channel, J. Fluid Mech., 135, 95–110. Rottman, J. W., J. E. Simpson, J. C. R. Hunt, and R. E. Britter (1985), Unsteady gravity current ﬂows over obstacles: Some observations and analysis related to the phase II trials, J. Hazardous Mater., 11, 325–340. Sequeiros, O. E. (2012), Estimating turbidity current conditions from channel morphology: A Froude number approach, J. Geophys. Res., 117, C04003, doi:10.1029/2011JC007201. Sequeiros, O. E., B. Spinewine, R. T. Beaubouef, T. Sun, M. H. Garcia, and G. Parker (2010a), Characteristics of velocity and excess density proﬁles of saline underﬂows and turbidity currents ﬂowing over a mobile bed, J. Hydraul. Eng., 136, 412–433. Sequeiros, O. E., B. Spinewine, R. T. Beaubouef, T. Sun, M. H. Garcia, and G. Parker (2010b), Bedload transport and bed resistance associated with density and turbidity currents, Sedimentology, 57, 1436–1490. Shin, J., S. Dalziel, and P. F. Linden (2004), Gravity currents produced by lock exchange, J. Fluid Mech., 521, 1–34. Stacey, M. W., and A. J. Bowen (1988), The vertical structure of turbidity currents and a necessary condition for self-maintenance, J. Geophys. Res., 94, 3543–3553. Tokyay, T., G. Constantinescu, and E. Meiburg (2011a), Lock exchange gravity currents with a high volume of release propagating over an array of obstacles, J. Fluid Mech., 672, 570–605. Tokyay, T., G. Constantinescu, E. Gonzales-Juez, and E. Meiburg (2011b), Gravity currents propagating over periodic arrays of blunt obstacles: Effect of the obstacle size, Journal of Fluids and Structures, 27, 798–806, doi:10.1016/j.jﬂuidstructs.2011.01.006. Tokyay, T., G. Constantinescu, and E. Meiburg (2012), Tail structure and bed friction velocity distribution of gravity currents propagating over an array of obstacles, J. Fluid Mech., 694, 252–291. Wynn, R. B., N. H. Kenyon, D. G. Masson, D. A. V. Stow, and P. P. E. Weaver (2002), Characterization and recognition of deep-water channellobe transition zones, AAPG Bull., 86, 1441–1462. Xu, J. P., M. A. Noble, and L. K. Rosenfeld (2004), In-situ measurements of velocity structure within turbidity currents, Geophys. Res. Lett., 31, L09311, doi:10.1029/2004GL019718. Yuksel, O. A., G. Constantinescu, and T. Tokyay (2012), Propagation of a gravity current in an aquatic canopy: Insights from large Eddy simulations, in River Flow 2012 International Conference on Fluvial Hydraulics, San Jose, Costa Rica, edited by R. E. Murillo Munoz, pp. 295– 303, CRC Press, Taylor and Francis Group, N. Y.

TOKYAY ET AL.

C 2014. American Geophysical Union. All Rights Reserved. V

2768

Lock-exchange gravity currents with a high volume of release ...