J. Fluid Mech. (2012), vol. 703, pp. 279–314. doi:10.1017/jfm.2012.213

c Cambridge University Press 2012

279

Internal bores: an improved model via a detailed analysis of the energy budget Zachary Borden1 , Eckart Meiburg1 † and George Constantinescu2 1 Department of Mechanical Engineering, University of California at Santa Barbara,

Santa Barbara, CA 93106, USA 2 Department of Civil and Environmental Engineering, University of Iowa, Iowa City, IA 52242, USA

(Received 14 April 2011; revised 16 February 2012; accepted 30 April 2012; first published online 13 June 2012)

Internal bores, or internal hydraulic jumps, arise in many atmospheric and oceanographic phenomena. The classic single-layer hydraulic jump model accurately predicts the bore height and propagation velocity when the difference between the densities of the expanding and contracting layers is large (i.e. water and air), but fails in the Boussinesq limit. A two-layer model, which conserves mass separately in each layer and momentum globally is more accurate in the Boussinesq limit, but it requires for closure an assumption about the loss of energy across a bore. It is widely believed that bounds on the bore speed can be found by restricting the energy loss entirely to one of the two layers, but under some circumstances, both bounds overpredict the propagation speed. A front velocity slower than both bounds implies that, somehow, the expanding layer is gaining energy. We directly examine the flux of energy within internal bores using two- and three-dimensional direct numerical simulations and find that although there is a global loss of energy across a bore, a transfer of energy from the contracting to the expanding layer causes a net energy gain in the expanding layer. The energy transfer is largely the result of turbulent mixing at the interface. Within the parameter regime investigated, the effect of mixing is much larger than non-hydrostatic and viscous effects, both of which are neglected in the two-layer analytical models. Based on our results, we propose an improved two-layer model that provides an accurate propagation velocity as a function of the geometrical parameters, the Reynolds number, and the Schmidt number. Key words: atmospheric flows, geophysical and geological flows, gravity currents

1. Introduction A bore, or a hydraulic jump, is an open channel flow phenomenon in which a sharp rise in height propagates along a fluid’s surface. Tidal bores, where the leading edge of the tide travels up a river against the current, are the best-known naturally occurring examples. An excellent description of the world’s major tidal bores can be found in Simpson (1997). The propagation of tidal bores is well-described analytically by single-layer hydraulic theory. Conserving mass and balancing the change in momentum flux with the change in hydrostatic pressure in a control volume surrounding a bore, as shown in figure 1(a), reveals a solution for the conditions † Email address for correspondence: [email protected]

280

Z. Borden, E. Meiburg and G. Constantinescu (b) p u

(a)

U2 pd U

hf U

ha

H

U1

hf U

U1

ha

F IGURE 1. Idealized geometry of: (a) a one-layer bore; and (b) a two-layer bore in a reference frame moving with the bore.

downstream of the bore (U1 and hf ) as a function of the upstream variables (U and ha ). This solution requires a loss of energy across the bore (Landau & Lifshitz 1959, p. 413). Bores can propagate along an interface between any two fluids, not just water and air as with a tidal bore. If a bore forms at the interface between two fluids of similar densities, such as fresh water and salt water, it is called an internal bore. The most studied naturally occurring internal bore is the ‘Morning Glory’ cloud, which occurs off of the northern coast of Australia (Clarke, Smith & Reid 1981). Bores similar to the Morning Glory have also been observed off the eastern coast of Florida, in the USA, and are created by the interaction of sea-breeze fronts and temperature inversion layers (Wakimoto & Kingsmill 1995). Likewise, thunderstorm outflows can also perturb an inversion layer and create atmospheric bores (Haase & Smith 1984). In marine environments, internal bores can arise from the breaking of internal waves (Leichter et al. 1996), by the interaction of the tides with ocean floor topography (Morozov et al. 2002), and by gravity current flows past submarine obstacles (Gonzalez-Juez, Meiburg & Constantinescu 2009; Gonzalez-Juez & Meiburg 2009; Gonzalez-Juez et al. 2010). These marine bores can have a substantial impact on the local environment. In lakes, internal bores have been proposed as a mechanism by which energy is transferred from basin-scale waves to solitons of a much shorter wavelength (Horn, Imberger & Ivey 2001), and in the deep ocean, Hosegood & van Haren (2004) argue that internal bores are a dominant mechanism for sediment resuspension and transport up the continental slope. In the first analytical models of internal bores, single-layer hydraulic theory was applied to a shallow layer beneath a semi-infinite fluid of finite density. This approach produced the same expression for front velocity as the single-layer hydraulic jump model, but with gravity replaced by a reduced gravity term to account for the reduced role of buoyancy (Turner 1973, p. 67). As the height of the fluid layer ahead of a bore, ha , goes to zero, the bore should behave as a gravity current with a welldefined propagation speed given by Benjamin (1968). The early single-layer theory, however, predicted that a bore’s velocity becomes infinite in this limit. Problems in this limit frequently arise in the study of thin, viscous films running down a slope. Analytical models describing the film’s propagation require a precursor film to capture an accurate front velocity (Zhou et al. 2005). Yih & Guha (1955) argue that when the less dense layer has a finite density, it can distribute momentum and energy across a bore. A two-layer model, which incorporates the upper, less dense fluid layer, should therefore be used. In their two-layer model, the mass in each layer is conserved in a control volume surrounding the bore, and the total loss of momentum across the bore is balanced by the change in hydrostatic pressure (figure 1b).

Internal bores: an improved model via a detailed analysis of the energy budget

281

The two-layer model introduces three new variables: the velocity in the upper layer downstream of the bore, U2 , the channel height, H, and the dynamic pressure jump across the bore, 1p = pu − pd . These additional variables require a fourth jump condition for a solution, one which relates the upstream and downstream energy fluxes. Finding the correct energy jump condition has been the focus of several previous investigations, which are discussed below, but remains an open question. In the present investigation, we directly compute the kinetic energy of internal bores through two- and three-dimensional direct numerical simulations in order to gain physical insight into the correct energy jump condition. A similar approach has been successfully applied to study the exchange of energy in gravity currents (Ungarish & Huppert 2006, 2008), and to evaluate gravity current front conditions (Ungarish 2008, 2009). First, we discuss existing analytical models for internal bores and proposed energy jump conditions in § 2. We then present the governing equations and numerical methods used for our simulations in §§ 3 and 4. Section 5 presents a comparison of our simulated bores with experimental results and analytical models. Subsequently, we employ the simulations in order to understand the energy balance across internal bores in § 6. In particular, we demonstrate that the lower, more dense layer usually gains energy as it expands across the internal bore. Finally, we propose a novel jump condition for the energy across the bore in § 8, and compare its predictions with those of existing models. 2. Existing theoretical models An analytical description of internal bores must bridge the gap between smallamplitude gravity waves and gravity currents. As the change in the height of the dense fluid layer across a bore approaches zero and the channel height becomes infinite, the bore should behave as an internal wave and propagate with the velocity U = (g0 ha )1/2 , where g0 = g(ρ1 − ρ2 )/ρ1 is the reduced gravity and ha is the thickness of the lower layer ahead of the bore. As the difference in layer height becomes large, the bore’s front velocity should approach the speed of a gravity current proposed by Benjamin (1968):

H − hf 2H − hf U
, (2.1)
1/2 = hf hf + H g0 hf

where hf is the height of the gravity current. Early analytical models of internal bores treat bores as single-layer hydraulic jumps. A control volume is constructed around the bore (figure 1a), and mass and momentum is conserved, which is written as Uha = U1 hf ,

U 2 ha + 21 g0 h2a = U12 hf + 12 g0 h2f .

(2.2) (2.3)

The solution of two these equations leads to a relation between the bore’s front velocity, uf , and the bore’s height, hf /ha , which is  1/2  1 hf hf uf = +1 , (2.4) 2 ha ha where uf = U/(g0 ha )1/2 is a non-dimensional front velocity. Although (2.4) correctly reduces to the velocity of a gravity wave when hf = ha , it blows up as ha → 0 and is therefore inconsistent with the gravity current limit described above.

282

Z. Borden, E. Meiburg and G. Constantinescu

As an improvement, Yih & Guha (1955) proposed a two-layer model which included the mass and momentum transport in the upper, less dense layer. Using the control volume in figure 1(b), they argued that if the two layers were immiscible, the bore’s propagation could be described by conserving mass in each layer and momentum globally. In a reference frame moving with the bore, these expressions are written as U1 hf = Uha ,

Z 0

H

U2 (H − hf ) = U(H − ha ), Z H

2 pu + ρg(H − z) + ρU dz = pd + ρg(H − z) + ρu2d dz,

(2.5) (2.6) (2.7)

0

where U, U1 and U2 are averaged velocities in each layer, ρ1 and ρ2 indicate averaged densities, pu and pd represent the pressures at the top of the channel upstream and downstream of the bore, and ud denotes the local velocity downstream of the bore. Combining (2.5)–(2.7) in the Boussinesq limit, however, reveals that more information is needed to find the bore’s speed in the two-layer model:  

H 1 ha H − ha hf − ha U 2 + g0 h2f − h2a + 1pH = 0. − + (2.8) hf H − hf 2 ρ1 An additional expression is needed to relate 1pH , the pressure drop along the top of the channel, to the upstream and downstream conditions. Chu & Baddour (1977) and Wood & Simpson (1984) proposed that this fourth expression should be a jump condition relating the up- and downstream energy fluxes. Based on observed mixing in dye streak experiments, they argued that energy is conserved across a bore in the contracting layer. They applied Bernoulli’s equation with no head loss along a streamline in the upper layer, which produces " # 1 (H − ha )2 2 1pH = ρ1 1 − (2.9)
2 U 2 H − hf as the expression for the pressure drop across a bore. Combining (2.8) and (2.9) gives an expression for the front’s velocity: (

2 )1/2 g0 hf hf + ha H − hf
Uws = , (2.10) H h2f − 3ha hf + 2ha H

which we refer to hereafter as the WS model. Klemp, Rotunno & Skamarock (1997) derived a model with a better fit with available experimental data and gravity current theory by using the alternative jump condition, a conservation of energy in the lower layer. The authors applied Bernoulli’s equation along a streamline in that layer, and obtained  
1 h2a 1pH = ρˆ 1 − 2 U 2 − g0 hf − ha , (2.11) 2 hf (

)1/2 g0 h2f 2H − ha − hf H − hf 
 Ukrs = (2.12) H h2f + H hf + ha − 3hf ha

Internal bores: an improved model via a detailed analysis of the energy budget

283

as expressions for the pressure drop and front velocity, which we hereafter refer to as the KRS model. For the purpose of subsequent discussion and comparison with simulation results, we express the WS and KRS models in dimensionless form as 1/2  R (1 + R) (1 − Rr)2 , (2.13) uws = R2 r − 3Rr + 2  2 1/2 R [2 − r (1 + R)] (1 − Rr) ukrs = , (2.14) R2 r − 3Rr + R + 1 where u = U/ (g0 ha )1/2 , R = hf /ha and r = ha /H. We can now examine these models’ behaviour in the internal wave and gravity current limits. The internal wave limit is recovered when r → 0 and R = 1. In this limit, both models predict the correct wave speed U = (g0 ha )1/2 . Revealing the gravity current limit, however, requires that we change the way we non-dimensionalize (2.13) and (2.14). By scaling the propagation velocity with (g0 hf )1/2 and introducing a new dimensionless variable α = hf /H, (2.13) and (2.14) become  1/2 Uws (α + r) (1 − α) , (2.15)
1/2 = α 2 − 3αr + 2r g0 hf   Ukrs α (2 − α − r) (1 − α) 1/2 . (2.16)
1/2 = α 2 − 3α + α + r g0 hf In the gravity current limit, r → 0, which simplifies (2.15) and (2.16) to   1 − α 1/2 Uws ,
1/2 = α g0 hf  1/2 Ukrs (2 − α) (1 − α) = .
1/2 α+1 g0 hf

(2.17) (2.18)

The propagation velocity given by Ukrs in (2.18) is the same equation derived by Benjamin (1968) for the front velocity of a gravity current, but the velocity in (2.17) grows unbounded as the channel depth goes to infinity (α → 0). Shortly after the publication of the KRS model, Li & Cummins (1998) developed a more general relation for the energy lost across an internal bore, a relation which encompassed both the WS and KRS models. Instead of assuming that the energy lost across an internal bore is dissipated entirely in one layer, they provided relationships between the energy lost in either the upper or the lower layer, and the front velocity. For a given amount of energy loss in one of the layers, the bore’s velocity can be either of   (R − 1)(1 − r) R2 r − 3Rr + 2 2 e˙ u = − uf − uf + 1 + R , (2.19) 2 R (1 − Rr)2  2  (R − 1) R r − 3Rr + R + 1 2 e˙ l = − uf uf + Rr + r − 2 , (2.20) 2 R2 (1 − Rr)

0 3/2 5/2 ˙ where e˙ = E/(ρg ha ) is the non-dimensional rate of change of energy per unit spanwise length. We refer to these equations as the LC relation. From this relation, we

284

Z. Borden, E. Meiburg and G. Constantinescu 0.4 0.2 0 –0.2 –0.4 –0.6 1.1

1.2

1.3

1.4

1.5

1.6

F IGURE 2. Overall rate of energy loss and rate of change of energy in each layer given by the LC relation, shown as a function of the front velocity for R = 2 and r = 0.1. The WS and KRS front velocities are recovered when e˙ u or e˙ l is zero (from Li & Cummins 1998). WS

3.0

KRS 2.5

2.0

1.5

1.0

1

2

3

4

5

6

7

F IGURE 3. Non-dimensional front velocity plotted as a function of R = hf /ha for r = 0.02. Here, the dashed line represents the front velocity predicted by WS and the solid line corresponds to KRS. The squares denote experimental data obtained by Wood & Simpson (1984) at many Reynolds numbers from Re = 500 to Re = 3800.

can recover the WS model by setting e˙ u = 0, and the KRS model by setting e˙ l = 0 (figure 2). Assuming that neither of the layers can gain energy, the WS and KRS models represent the upper and lower bounds on the front velocity of an internal bore. Comparing the predictions of the WS and KRS models with experimental data from Wood & Simpson (1984) (figure 3) shows that the KRS model is the better fit, but it slightly overpredicts the bore’s velocity. For front velocities smaller than ukrs , the LC relation predicts a positive e˙ l – an increase of energy in the lower layer across the bore (figure 2). The existence of an energy gain in the lower layer is controversial. Baines (1995) and Wood & Simpson (1984) assumed that a gain of energy in either of the layers

Internal bores: an improved model via a detailed analysis of the energy budget

H

285

uf h0

z

hf

x

ha

F IGURE 4. Schematic of a dam-break problem in a laboratory reference frame. The dashed line is the initial condition and the solid line is the state after some time has passed.

was unphysical and argued against Yih & Guha’s original internal bore jump condition, which has not been discussed here, because it implied a slight energy gain in the upper layer. Klemp et al. (1997) noticed evidence for a small energy gain in the lower layer in some of their modelling studies and argued that shear at the interface between the two layers could transfer enough energy to the lower layer to overcome dissipation. Their studies, however, were conducted on a coarse numerical mesh (300 × 100) and used a turbulence closure scheme to represent viscous energy loss, leaving them unable to calculate the energy gain accurately (Klemp, Rotunno & Skamarock 1994). Finally, Li & Cummins (1998) remark that a gain of energy in the lower layer cannot be supported because of a lack of experimental data. Further study is needed. 3. Governing equations To simulate an internal bore, we model a dam break, where a finite reservoir of dense fluid is suddenly released into a two-layer channel (figure 4). To simplify our modelling, we invoke the Boussinesq approximation. Earlier experimental work on gravity currents by Huppert & Simpson (1980) found that gravity currents produced by a dam break reach a steady front velocity after an initial transient phase. Internal bores behave similarly, which we show in detail in § 4.2 by investigating the temporal evolution of the bore’s height and front velocity. In a two-dimensional simulation, which forms the basis for most of our parametric studies, the velocity and pressure are governed by the Navier–Stokes equations, and the density evolves by an advection–diffusion equation. These are written as

∇ · u = 0,

(3.1)

∂u 1 2 + u · ∇u = −∇p + ∇ u + ρeg , (3.2) ∂t Re ∂ρ 1 + u · ∇ρ = ∇ 2 ρ, (3.3) ∂t Re Sc where u = (u, w)T denotes the velocity vector, p the pressure, ρ the density, and eg the unit vector in the direction of gravity. Because we consider only horizontal channels, eg = (0, −1)T . All terms in (3.1)–(3.3) have been made dimensionless by the undisturbed fluid height h˜ a , the upper-layer density ρ˜2 , and the buoyancy velocity u˜ b , which we define as  1/2 u˜ b = g˜ 0 h˜ a . (3.4)

286

Z. Borden, E. Meiburg and G. Constantinescu

Here, a tilde denotes a dimensional quantity. The non-dimensional pressure p and density ρ are given by p=

p˜ , ρ˜2 u˜ 2b

ρ=

ρ˜ − ρ˜2 . ρ˜1 − ρ˜2

(3.5)

The two non-dimensional parameters appearing in (3.2) and (3.3) are the Reynolds number, Re, and the Schmidt number, Sc, which we define as u˜ b h˜ a ν˜ , Sc = , (3.6) ˜ ν˜ D ˜ is the molecular diffusivity. where ν˜ is the kinematic viscosity and D Our two-dimensional simulations will unavoidably neglect some three-dimensional effects, especially in bores that are turbulent. For example, viscous energy dissipation will be stronger in three-dimensional bores because of additional instabilities in the spanwise direction. We therefore conduct some fully three-dimensional simulations to verify that our two-dimensional results can be generalized to three dimensions. Our three-dimensional equations, methods, and results are presented in § 7. Re =

4. Numerical methods For the two-dimensional dam-break simulations, we use the streamfunction-vorticity code developed by H¨artel, Meiburg & Necker (2000). Those authors extensively validated the code by applying it to test problems with known solutions and obtained excellent agreement. We additionally verified that the code conserves energy by tracking the volume integral of the kinetic energy, the volume integral of the potential energy, and the time integral of the spatially integrated dissipation field. After t = 40 non-dimensional time units, the runtime of our longest simulation, the sum of these three terms is constant to within 0.5 % of the original potential energy in the computational domain. Our code uses Fourier sine or cosine expansions for streamwise discretization and sixth-order compact finite differences for normal spatial discretization. For boundary conditions, we typically apply a perfect-slip condition rather than the more realistic noslip condition. We choose this condition to more accurately approximate the theoretical control-volume models, which do not account for dissipation along the boundaries. To be thorough, we also run simulations with no-slip boundaries in § 6.5 to examine their influence. The slip walls on the domain’s streamwise boundaries can produce reflections, but we choose the domain’s streamwise length, 2L, to be sufficiently long such that our simulations terminate well before the front of the bore reaches the downstream boundary. The precise value of L varies depending on R and r, but is typically of the order of L ≈ 60. The domain’s height, H, is a direct function of r. Because we non-dimensionalized the governing equations on the length scale ha , the height is defined as H = 1/r. Our use of spectral methods necessitates a smooth initial condition. We use a twodimensional Gaussian function with a width of 0.07 to smooth our initial density field. For the fine grids we use, the width of the Gaussian function does not measurably affect the front velocity or height during the quasi-steady phase of a simulation. In their laboratory experiments, Wood & Simpson (1984) obtain a range of Reynolds numbers from Re = 500 to Re = 3800. We usually use a resolution of 8193 × 513 grid points, allowing us to simulate their entire range of Reynolds

Internal bores: an improved model via a detailed analysis of the energy budget

287

H L

F IGURE 5. Density field of a representative simulation with snapshots taken every two or four dimensionless units of time. Here, R = 2.43, r = 0.2, Re = 3500 and Sc = 1. The domain half-length, L = 30, is sufficiently long such that the bore never reaches the end of the domain during the simulation. The domain extends to x = −30 to the left, but is truncated in this figure.

numbers. We also performed some simulations at a much more refined resolution of 16 384 × 2048, which made our grid Reynolds number of order unity. But we found that our energy results, discussed below, were qualitatively identical, and quantitatively different by only a few per cent. Several density fields taken at different times in a representative simulation are shown in figure 5. To compare the simulated bore’s front velocity with theoretical models and experimental data, we need to find the speed of the bore uf as a function of r, R, and the energy fluxes in each layer. The following sections will describe how we calculate these quantities. 4.1. Finding the propagation velocity

To find a bore’s front velocity, we track the front position, xf (t), and differentiate it with respect to time. Tracking the front position would be simple with a well-defined shock, but the front of an internal bore is much smoother, and so its position is ambiguous. Wood & Simpson (1984) found that different front tracking methods yielded different front velocities. Consequently, we examine the two most popular methods. Our first method defines the front’s location as the first streamwise position, from the right, where the height of the ρ = 0.5 isopycnal is above a critical value, which we define as hc =

h0 − ha + ha . 4

(4.1)

In the other method, we mark the fluid initially behind the dam with ‘dye’. The dye concentration, cd , is governed by the advection–diffusion equation ∂cd 1 + u · ∇cd = ∇ 2 cd , ∂t Re Sc

(4.2)

288

Z. Borden, E. Meiburg and G. Constantinescu (b)

(a) 50

(i)

1.5 40 (i)

30 20

1.0 (ii)

(ii)

0.5

10 0

10

20

30

40

0

10

20

30

40

Time

Time

F IGURE 6. Comparison of the two methods for determining the front position (a) and velocity (b): (i) critical height; and (ii) dye method. The front velocity we extract from the critical height method is shown with a dotted line in (b). R = 2.3, r = 0.1, Re = 3500 and Sc = 1.

but is not coupled into (3.2). It therefore has no influence on the flow dynamics. We define the front’s location as the rightmost position of the cd = 0.5 dye concentration contour, where the initial dye concentration in the dense fluid is unity. For both methods, a second-order polynomial is used to interpolate the front’s location between gridlines. Despite interpolation, the position data are uneven as the front’s location moves between gridlines. We therefore have to smooth the data with local averaging before differentiating to avoid noisy results. In a typical simulation, the front velocity from both methods approaches an asymptotic value after an initial transient and becomes steady (figure 6). But the two velocities differ significantly. Wood & Simpson (1984) noticed this discrepancy in their experiments and attributed it to viscous effects at the bottom wall. Our simulations have no viscous effects at the bottom wall because of the free-slip boundary conditions, but we still observe that the bore moves into the clear fluid. The dye does not reliably track the front’s location. We therefore use the threshold height measurement, which is also consistent with the experimental studies of Baines (1984) and Wood & Simpson (1984). To calculate a single value of the front velocity, we take the average from t = 25 to t = 35. 4.2. Finding the height of the jump Because we are simulating a dam break to create an internal bore, the height of the bore, hf , and therefore R, cannot be specified a priori; it must be measured. For low Re and high Sc, the interface between the two layers is clearly defined, and the local height h(x) can be determined by inspection. When mixing and diffusion occur, however, the interface becomes diffuse and harder to characterize. We therefore define the local height as Z H h(x) = ρ (x, z) dz, (4.3) 0

which describes the interface as if the velocity field were suddenly frozen and all dense fluid allowed to settle. Usually in experimental studies of internal bores, hf is determined by inspection at one location (Baines 1984), but because our study is numerical, a few alternative procedures become available. We explore three options for calculating hf in order to

Internal bores: an improved model via a detailed analysis of the energy budget

289

(b) 2.5

(a) 2.0

2.0

h

1.5

1.5

Measured Flux Gate

1.0

1.0 0

10

20

30

40

0

20

40

60

x

Time

F IGURE 7. Comparison of the three methods for finding hf . (a) How the height calculated using each method varies over time for a simulation with R ≈ 2, r = 0.1, Re = 3500 and Sc = 1. The direct measurement and flux methods eventually converge to a nearly steady value. (b) The calculated heights overlaid on the local height at time t = 32. The gate and measured heights are indistinguishable.

identify the most appropriate, and refer to them as the direct measurement, the flux, and the gate methods. In the direct measurement method for finding hf , we take a spatial average of the local height between the bore’s front and the rarefaction wave propagating into the initial reservoir. This region expands over time, so we visually determine the edges of the region for five times and use a linear fit to interpolate intermediate values. In the flux method, we define a control volume encompassing the lower layer with the left edge pinned at x = 0, and the right edge attached to the tip of the bore. The bore’s height is then found by dividing the flux of lower-layer fluid into the control volume by the control volume’s streamwise rate of expansion. The flux of lower-layer fluid into the control volume at x = 0 is given by the expression Z H φ= (4.4) (ρu)x=0 dz, 0

and the amount entering from the right due to the control volume’s expansion is uf ha . By dividing the sum of these two fluxes by the control volume’s rate of streamwise expansion, the bore’s front velocity, we find that the average height is given by hf =

φ + uf ha . uf

(4.5)

Finally, in the gate method, we define hf as the local height of the lower layer at the dam location. We do not expect this method to produce a reliable height measurement because of temporal fluctuations, but include it for completeness. A comparison of the three methods is shown in figure 7. Each method produces slightly different results which, with the exception of the gate method, attain a quasi-steady value after t = 25. Both the direct measurement and the gate method’s results are the most accurate in a visual comparison with the local height profile in figure 7(b), but the gate method is prone to temporal oscillations because of its narrow spatial sampling. The flux method overpredicts the bore’s height because of undular waves at the bore’s front. These waves increase in both number and amplitude over time, acting as a sink of mass that is not accounted for in (4.5). Because of the

290

Z. Borden, E. Meiburg and G. Constantinescu (b)

(a) WS

1.40

KRS

1.35

1.40

WS

1.35

KRS

1.30

1.30

1.25

1.25

1.20

0

1000

2000

3000

4000

5000

1.20

0

1

2

3

4

F IGURE 8. Influence of the flow parameters Re (a) and Sc (b) on the propagation velocity for r = 0.1 and R ≈ 1.9. Here, the symbols indicate our simulation results, the solid line represents the KRS model, and the dashed line represents the WS model. In the Reynoldsnumber study, Re Sc = Pe = 3500 and in the Schmidt-number study, Re = 3500.

Ideal U2

Real U2

U2

hf U U1

U1

U1 U

F IGURE 9. Internal bore with an expanded view near the interface between the upper and lower layers. The ideal flow field assumed in the analytical models contains a sharp change in velocity at the interface. Realistically, the velocity must transition gradually over some distance δ.

disadvantages of the flux and gate methods, we use the direct measurement method in our analyses, and take an average height from t = 25 to t = 35 to further reduce temporal oscillations. 5. Parameter study Here, we present the results of several parameter studies that show the influence of Re, Sc, R and r on a bore’s propagation velocity, uf . Below, these studies will be used to compare our improved model with the existing two-layer models. We first examine the influence of the flow parameters Re and Sc. Figure 8(a) shows that above Re ≈ 1000, the front velocity is independent of the Reynolds number. But as the Reynolds number decreases, so does the front velocity. To understand this decrease in the front velocity, consider the shear layer between the upper and lower layer for a finite Reynolds number (figure 9). In low-Reynolds-number flow, the thickness δ of the shear layer increases with decreasing Reynolds number. Boundary layer theory shows that for laminar flow, δ ∼ Re−1/2 . As δ increases, the interface

Internal bores: an improved model via a detailed analysis of the energy budget (a)

291

(b) 1.6 1.6 1.4 1.4 1.2

1.2

1.0

1.0 1.2

1.4

1.6

1.8

2.0

R

2.2

2.4

2.6

0

0.05

0.10

0.15

0.20

0.25

r

F IGURE 10. Influence of the geometrical parameters, (a) R and (b) r, on the propagation velocity for Re = 3500, Sc = 1 and different values of r and R. Here, the filled symbols (, , H) represent our simulation results, the open circles (◦) represent experimental data from Wood & Simpson (1984), the solid lines represent the KRS model, and the dashed lines represent the WS model.

•

between the two layers of a bore becomes less well-defined, which means that the way mass and momentum conservation are stated in (2.5)–(2.7) becomes increasingly questionable. Klemp et al. (1997) derived a two-layer model in which the shear layer has a finite thickness and found that to first order, the front velocity decreased with increasing interface thickness, consistent with our observations. The Schmidt number’s influence on the front velocity of our simulated internal bores is shown in figure 8(b). Previous work by Necker et al. (2005) on gravity currents suggests that the large-scale flow properties of gravity currents, such as front velocity and current height, depend on Sc only when it is smaller than 1. Bonometti & Balachandar (2008) also find that the dynamics of gravity currents are
independent of the Schmidt number when Sc > 1, but only when Re > O 104 . For currents with smaller Re, the value beyond which Sc no longer affects global flow properties increases. We observe that a bore’s Schmidt-number dependence is similar to that of gravity currents. The front velocity is strongly dependent on the Schmidt number for low Sc but appears to level off at Sc ≈ 2. Simulations at higher Sc would be necessary to strengthen this claim, but numerical stability issues prevented us from completing these simulations in a reasonable amount of time. Consider that there are two boundary layers at the interface between the upper and lower layers: a momentum boundary layer and a density boundary layer. We suppose that the large-scale flow behaviour will be dominated by the thicker of these two layers. If Sc > O(1), the viscous boundary layer is thicker, and the magnitude of δ will be governed by the Reynolds number; the Schmidt number will not influence the front velocity. If the Schmidt number is small, however, the density boundary layer is thicker and δ will be determined primarily by Sc. Next, we examine the influence of the geometrical parameters, R and r, which were extensively studied by Wood & Simpson (1984), Klemp et al. (1997), and Li & Cummins (1998) when they developed their two-layer models. From our parameter study, we observe that as R increases, so does the propagation velocity (figure 10a). The increase in front velocity is not surprising because as the height of a bore increases, the horizontal pressure gradient driving its movement becomes

292

Z. Borden, E. Meiburg and G. Constantinescu 4 3

z 2 1 –10

–5

0

5

10

15

20

25

30

35

40

x

F IGURE 11. Height profile of bores with R = 1.33, 1.52, 1.69, 1.85, 2.04, 2.2, 2.38, 2.57, 2.67 in the reference frame of figure 1. Here, x = 0 corresponds to the front of the bore.

larger as a result of the increased hydrostatic pressure behind the front. But beyond some critical value of R, the front velocity begins to decrease with increasing R, which can be seen in the WS and KRS model predictions for r = 0.2. As the downstream thickness of an internal bore occupies a greater portion of the channel depth, the pressure gradient across the top of the channel required to accelerate the upper-layer fluid across the bore’s front increases. For a sufficiently large bore, increasing R will cause this opposing pressure gradient to grow more rapidly than the hydrostatic pressure gradient driving the bore’s propagation, and the front velocity will decrease. We also observe that the structure of the front changes with R. For small values of R below 1.9, the bore is undular and the front resembles a train of gravity waves. This can be seen most clearly with the smallest bore in figure 11. As R increases, turbulent mixing behind the first crest eliminates the waves. As r increases (as the channel height decreases), the dynamic pressure gradient along the top of the channel that opposes the bore’s propagation increases to accelerate the upper-layer fluid across the bore. This opposing pressure gradient causes a slower propagation (figure 10b). Similar to the previous experimental and numerical studies by Wood & Simpson (1984) and Cummins (1995), our simulated bores are not propagating as fast as the WS and KRS models predict they should. To discover the reason, we next examine a representative bore in more detail. 6. Energy transfer mechanisms In order to understand the mechanics of internal bores better, we examine a representative simulation with R = 1.85, r = 0.1, Re = 3500 and Sc = 1. All observations of this bore are made at time t = 32, well after the bore’s front velocity and height have reached a quasi-steady value. For our analysis, we define x = 0 to correspond to the bore’s front, and change the reference frame into one moving with the front of the bore to match the reference frame in the WS and KRS models. Figure 12 shows that in this new reference frame, the slip boundaries on the top and bottom walls produce nearly uniform and horizontal flow in each layer both upstream and downstream of the bore. The total pressure minus the hydrostatic component from the upper layer is shown in figure 13 for our representative bore. We are using a streamfunction–vorticity formulation of the Navier–Stokes equations in our simulations, which means that the pressure is never explicitly determined. To calculate the pressure field, we use the

Internal bores: an improved model via a detailed analysis of the energy budget

293

6 5 4

z 3 2 1 0 –15

–10

–5

0

5

10

15

20

25

30

35

40

x

F IGURE 12. Velocity field overlaid on the density profile shown in a reference frame moving with the bore. The velocity field is shown only at every 60th point in x and 20th point in z. Both ahead of and behind the bore, the velocity profile is horizontal and uniform in each layer. Across the bore, the flow accelerates in the upper layer and slows down in the lower layer. Here, x = 0 corresponds to the front of the bore. 6 5 4

z 3 2 1 0 –15

–10

–5

0

5

10

15

20

25

30

35

40

x

F IGURE 13. Pressure field and contours for the representative bore. Here, darker areas correspond to higher pressures. The thick black line indicates the local height of the bore.

approach described in Fletcher (1991) and solve the Poisson equation   ∂u ∂w ∂w ∂u ∂ρ ∇ 2p = 2 − − , ∂x ∂z ∂x ∂z ∂z

(6.1)

with the condition ∂p/∂x = 0 along the left and right boundaries, and ∂p/∂z = −ρ along the top and bottom boundaries. Along the bottom of the channel, the pressure increases as we move from left to right across the bore because of the increased hydrostatic pressure force of the more dense, lower-layer fluid occupying a greater fraction of the channel depth. This groundlevel pressure jump was also recorded by Clarke et al. (1981) when an atmospheric bore passed overhead. The low-pressure regions atop the undular crests near the bore’s front are the result of Bernoulli’s theorem; the increased velocity in that region causes a drop in pressure.

294

Z. Borden, E. Meiburg and G. Constantinescu

In the atmosphere, the low pressure would cause a corresponding drop in temperature, which could cause the formation of clouds, depending on the humidity. These regions of low pressure are therefore consistent with the long roll clouds sometimes observed in the Morning Glory atmospheric bore. As we stated in the introduction, our goal is to investigate the energy jump conditions needed for closure in the two-layer models by tracking the evolution of the kinetic energy in our simulated internal bores. The local time rate of change of kinetic energy is governed by
∂E ∂ + u · ∇E = ρeg · u − u · ∇p + 2Re−1 ui Sij + Φ , ∂t ∂xj |{z} | {z } | {z } | {z } | {z } |{z} (a)

(b)

(c)

(d)

(e)

(6.2)

(f )

where E is the non-dimensional kinetic energy defined as ui ui /2, Sij denotes the rate-of-strain tensor, and Φ represents the viscous dissipation (Panton 2005, p. 89). From left to right, the terms in (6.2) represent: (a) the local time rate of change of kinetic energy; (b) the convective flux of kinetic energy; (c) the work performed by body forces; (d) the work performed by pressure forces; (e) the viscous diffusion of momentum (equivalently, the work performed by viscous stresses); (f ) the viscous loss of kinetic energy to heat. In our two-dimensional system, the viscous diffusion of momentum is given by
∂ 2Re−1 ui Sij ∂xj # "    2   ∂w ∂u ∂w 2 ∂u 2 −1 2 2 +2 + + u∇ u + w∇ w , = Re 2 + ∂x ∂z ∂z ∂x and the viscous dissipation, Φ, is given by "    2  # ∂u 2 ∂w ∂u ∂w 2 −1 −1 Φ = −2Re Sij Sij = −Re 2 +2 + + . ∂x ∂z ∂z ∂x

(6.3)

(6.4)

To directly compare our simulations with the existing two-layer models, we integrate the terms in (6.2) in a control volume extending across the bore and encompassing either the whole channel or each layer separately. We choose a control volume that extends from the edge of the domain ahead of the bore to some location xr downstream of the bore’s front. For a control volume spanning the entire channel height, we evaluate Z xr Z H At = a dz dx, (6.5) x0

0

where a is the per-unit-volume quantity we wish to integrate and x0 is the location of the wall ahead of the bore in the reference frame centred on the bore’s front. To

Internal bores: an improved model via a detailed analysis of the energy budget

295

isolate a single layer, we use the local height h(x) to define the limits of integration Z xr Z h(x) Al = a dz dx, (6.6) x0

Au =

Z

xr

x0

0

Z

H

a dz dx.

(6.7)

h(x)

Because the front velocity of the bore is slower than the velocity predicted by both the WS and KRS models, we expect to find a gain of energy in the lower layer across the bore, as well as a global loss of energy. We now examine how such a change in energy would appear in (6.2). The mechanical energy evolution equation states that pressure, body, and viscous forces can affect either the time rate of change of kinetic energy within the control volume (term a), or the flux of kinetic energy across the control-volume boundaries (term b). Because the front velocity and height of our bore are constant at the time of our analysis, we expect the time rate of change of kinetic energy (term a) to equal zero in any control volume we choose. In reality, this term will not exactly equal zero because of transient fluctuations in the mixing region between the two layers, but it is much smaller in magnitude than the convective flux of kinetic energy. In a control volume encompassing the entire channel height and extending from ahead of the bore to xr = 40 units downstream of the front, the convective flux of kinetic energy is 38 times greater in magnitude than the time rate of change of kinetic energy for the representative bore. In the momentum balance adopted in the two-layer bore models, (2.7), the change in momentum flux across a bore is caused by the horizontal pressure gradient that results from the change in hydrostatic pressure across the bore, as well as the change in pressure across the top of the channel. Taking the dot product of this momentum balance with the velocity field gives us the expected energy balance in the theoretical model, namely, that the change in the convective flux of kinetic energy is only the result of work done by the horizontal pressure gradient. If there were no energy loss across a bore in one of the layers, this expected energy balance would be satisfied exactly in that layer. In the framework of (6.2), term (b) integrated in a control volume spanning a single layer would exactly cancel the u ∂p/∂x component of term (d) integrated in the same control volume. When we say a layer is gaining or losing energy in the sense of WS and KRS, we really refer to any change in convective flux of kinetic energy that is not the result of the opposing horizontal pressure gradient. Within each layer, the flux of kinetic energy changes across the bore. For example, the opposing horizontal pressure gradient will decelerate fluid in the lower layer, and the layer’s convective flux of kinetic energy should decrease as well. But, if the convective flux decreases by more than the work by the horizontal pressure gradient dictates it should, we say the layer has ‘lost energy’. If the opposite is true, it has ‘gained energy’. It is important to note that this change in energy may or may not have anything to do with the actual viscous dissipation of energy into heat. Viewing the energy loss in this manner, we explore three mechanisms by which a layer can gain or lose energy. First, kinetic energy can diffuse and be dissipated because of viscous effects. Viscous mechanisms are irreversible and represent the only way that the total flux of energy, kinetic plus potential, can decrease across a bore. Second, non-hydrostatic effects, which are ignored by the two-layer control-volume models, could somehow convert kinetic energy to potential energy, or vice versa. The conversion cannot affect the total amount of energy, but if it affects the bore’s velocity,

296

Z. Borden, E. Meiburg and G. Constantinescu 0.01

0

–0.01

–0.02

Total Lower Upper

–0.03 –10

0

10

20

30

F IGURE 14. Change in the flux of kinetic energy due to viscous diffusion and dissipation plotted as a function of the position of the right-hand boundary of a control volume spanning both layers, or the upper or lower layer separately. The viscous diffusion of momentum does not affect the energy flux across a global control volume because of the slip boundaries, so the solid line represents the global loss of energy due to viscous dissipation.

it will appear as an energy loss or gain in the LC relation. Finally, mixing at the interface of the two fluid layers due to turbulence can also convert kinetic energy to potential energy by lifting denser fluid parcels above lighter ones. Just as with the non-hydrostatic effects, this would appear as a gain or loss of energy in the LC relation. We now discuss each of these three mechanisms in detail. 6.1. Viscous effects Both the WS and KRS models predict a global loss of energy across an internal bore. Because we are conducting direct numerical simulations, the only ‘real’ mechanism for energy loss is the viscous dissipation of kinetic energy into heat, which is represented by term (f ) in (6.2). In a global control volume extending to xr = 35, the total change in energy flux due to viscous dissipation in the representative simulation is e˙ t = −1.95 × 10−2 . The amount of dissipation is split fairly evenly between the two layers, with 48 % occurring in the upper layer and 52 % occurring in the lower. The other viscous term in (6.2), term (e), represents the viscous diffusion of momentum due to shear, that is to say, the work performed by viscous forces. Because we employ slip boundaries in our simulations, there is no momentum flux through the top and bottom walls. This term therefore cannot affect the total energy in a global control volume, it can only redistribute energy within it. The velocity gradient at the interface between the upper and lower layers will cause momentum to diffuse from the faster moving, upper-layer fluid into the slower moving lower layer. When we plot the sum of the two viscous terms as a function of control-volume size (figure 14), we find an increase of energy in the lower layer. The viscous diffusion of momentum transferred more energy into the lower layer than was lost by viscous dissipation in that layer. We also see in figure 14 that if the right-hand side of the control volume is located beyond the undular front of the bore (after xr = 15), the magnitude of the viscous terms grows approximately linearly with the size of the control volume. The WS

Internal bores: an improved model via a detailed analysis of the energy budget

297

0.12 Total Lower Upper

0.10 0.08 0.06 0.04 0.02 0 –10

0

10

20

30

F IGURE 15. Change in the flux of kinetic energy due to the non-hydrostatic terms plotted as a function of the position of the right-hand boundary of a control volume spanning both layers, or the upper or lower layer separately.

and KRS models’ predictions about the total energy losses, however, are independent of the control-volume size. The only way we can make comparisons between the simulations and the theoretical models is to pick a standard control-volume size, such as xr = 35, and apply it to all of our simulations. Fixing the size of the control volume will ensure that trends within the results of our simulations are accurate, but will leave a comparison with the analytical models undetermined to a constant. 6.2. Non-hydrostatic effects In internal bores, as with many other atmospheric and oceanographic currents, the streamwise length scales and velocities are much larger than the vertical length scales and velocities. Both the WS and KRS models therefore invoke the hydrostatic assumption, which states that the vertical pressure gradient is purely the result of the gravitational body force. If our simulated internal bores were perfectly hydrostatic, term (c) and the w ∂p/∂z component of term (d) in (6.2) would exactly cancel each other over the entire domain. The bore cannot be perfectly hydrostatic, however, because of the vertical acceleration of fluid at the front of the bore, combined with the presence of the low-pressure regions atop the undular crests of the bore (figure 13). It could be that the non-hydrostatic terms in (6.2), which are neglected in the WS and KRS models, are affecting the convective flux of kinetic energy. Any change in the convective flux that is not the result of the horizontal pressure gradient will appear as a gain or loss of energy in the framework of the theoretical models. Figure 15 shows the change in the convective flux of kinetic energy as a result of the imbalance between term (c) and the vertical component of term (d). We see that there are large fluctuations in a global control volume as its right-hand edge passes through the front of the bore, but downstream, the fluctuations die out and reveal only a slight energy gain. Furthermore, this energy gain is confined almost entirely to the upper layer. There is almost no change in the energy flux in the lower layer. Comparing the amount of energy change from the viscous and non-hydrostatic terms is difficult. For the representative control volume we have chosen, with xr = 35, the global change in the convective flux of kinetic energy due to the non-hydrostatic terms

298

Z. Borden, E. Meiburg and G. Constantinescu A

B

A

B

H

0

H

1

0

1

F IGURE 16. The effect of turbulent mixing on the density stratification downstream of an internal bore.

(˙et = 1.94 × 10−2 ) is almost identical in magnitude to the change in the convective flux of kinetic energy due to viscous dissipation. We could have chosen a different sized control volume though, say xr = 70, and the amount of viscous dissipation would be roughly doubled, while the energy change due to the non-hydrostatic terms would remain constant. The magnitude of the energy change resulting from non-hydrostatic effects also turns out to be much smaller than the effect of turbulent mixing, which we discuss next. We believe it is appropriate to conclude that non-hydrostatic effects can be ignored in the global energy balance. 6.3. Turbulent mixing

The final mechanism we explore which can change the convective flux of kinetic energy is turbulent mixing at the interface between the upper and lower layers. In the analytical control-volume models discussed earlier, it is assumed that the interface between the upper and lower layers is sharp and well-defined downstream of the bore’s front. We observe this sharp interface in our simulations very close to the bore’s front, but as we move farther downstream, interfacial shear results in turbulent mixing, which redistributes some lower-layer fluid into the upper layer, and vice versa. In the case of a shear layer in a homogeneous medium, the shear layer’s thickness, δ, grows linearly with downstream distance (Winant & Browand 1974). In the stratified case, however, Koop & Browand (1979) observed that stratified shear layers reach a steady-state thickness a certain distance downstream of their origin; an energy constraint puts bounds on the downstream growth. Consider the bore, and associated shear layer, in figure 16. At cross-section A–A, near the front of the bore, relatively little mixing has occurred and the layers’ interface is sharp. At cross-section B–B, much more mixing has occurred, resulting in a different density profile. The potential energy at A–A is given as PEA = h2f /2, whereas at B–B, if we assume a linear shear layer, the potential energy is given as PEB = h2f /2 + δ 2 /24. The potential energy at the mixed cross-section will always be greater than at the unmixed one, regardless of the value of δ. If the flux of potential energy increases with downstream distance, the flux of kinetic energy must decrease. Turner (1986) was able to derive an expression for the shear layer’s asymptotic thickness by conserving mass and momentum, equating the changes in potential and kinetic energy, and assuming a linear stratification in the

Internal bores: an improved model via a detailed analysis of the energy budget

299

0.2

0.1

0

–0.1 10 –4

10 –3

10 –2

10 –1

10 0

F IGURE 17. Calculated change in kinetic energy flux due to turbulent mixing in the lower (H), upper (N), and both ( ) layers as a function of the threshold density ρt .

•

velocity and density fields over the same distance. It is 2ρ2 (1U)2 . (6.8) g1ρ In (6.2), the equation governing the transport of mechanical energy, the effect of mixing is contained within the convective flux and horizontal pressure gradient terms. To extract it, we calculate an idealized flux of kinetic energy in the absence of mixing, and then subtract that quantity from the actual flux of kinetic energy. To calculate the idealized flux, we first determine the average velocity in each layer in the absence of mixing. We define a threshold density ρt and then find the average streamwise velocity for all x in cells where ρ > (1 − ρt ) for the lower layer, and ρ < ρt for the upper layer. The flux of kinetic energy is then computed with the equations δmax =

fl = 12 u¯ l (x)3 h(x),   1 3 1 fu = 2 u¯ u (x) − h(x) , r

(6.9) (6.10)

where u¯ l is the average velocity in the lower layer and u¯ u is the average velocity in the upper layer. Any difference between (6.9)–(6.10) and the actual flux of kinetic energy in each layer which is not the result of viscous or non-hydrostatic effects must be due to mixing at the interface. In this method, though, the magnitude of mixing is affected by the parameter ρt . We observe in figure 17 that the amount of mixing becomes independent of the threshold density once the threshold density is sufficiently small. We therefore always choose ρt = 2.5 × 10−3 . Turbulent mixing has a substantial impact on the flow. Figure 18 shows that there is an apparent global loss of energy due to mixing that is an order of magnitude larger than either the viscous or non-hydrostatic effects. The lower layer, however, gains energy because of turbulent mixing. For the representative simulation, this energy gain in the lower layer is six times greater than the energy gain due to viscous effects, and 40 times greater than the increase in energy due to non-hydrostatic effects. The change in kinetic energy flux due to turbulent mixing becomes independent of the control-volume size as long as the right-hand edge of the control volume is

300

Z. Borden, E. Meiburg and G. Constantinescu 0.2 0.1 0 –0.1 –0.2 Total Lower Upper

–0.3 –0.4 –0.5 –10

0

10

20

30

F IGURE 18. Change in the flux of kinetic energy due to turbulent mixing plotted as a function of the position of the right-hand boundary of a control volume spanning both layers, or the upper or lower layer separately.

sufficiently far behind the front of the bore, beyond xr = 25 for the present case. Koop & Browand (1979) noted that in stratified shear layers, the position of maximum thickness varies as a function of the Richardson number, which Koop & Browand (1979) defined as Ri = g1ρhi /ρa (1U)2 , where hi is the thickness of the shear layer 1 cm behind its front, and ρa is the average density. Their choice of length scale makes it difficult for us to compute a Richardson number because the value of hi is sensitive to the position where it is calculated. If we instead define hi = 0.1, which is the thickness of the density stratification upstream of the bore, Ri = 0.145 for our simulation. With this value of Ri, the position of maximum thickness of the shear layer should be located around xr = 30, which compares well with our finding that no more energy transfer due to mixing occurs after xr = 25. 6.4. Energy parameter study We have shown that turbulent mixing at the horizontal interface is responsible for the energy gain in the lower layer, and now we examine how the change in energy flux due to mixing is influenced by Re, Sc, R and r. The influence of the Reynolds number on the change in energy flux from mixing is shown in figure 19(a). For all values of Re, mixing causes a gain of energy in the lower layer, but a global loss of energy, which is much greater in magnitude than the actual energy loss due to viscous dissipation for Re > 1000. For smaller values of Re, viscous effects become important and are of the same order as those of turbulent mixing. For Re > 1000, the energy change in both layers due to mixing is a very weak function of the Reynolds number. We therefore expect, via the LC relation, that the front velocity should also be almost independent of the Reynolds number, which we earlier noted in figure 8. The results of the Schmidt-number study are also straightforward to interpret (figure 19b). As Sc decreases below O(1), molecular diffusion becomes the dominant mechanism controlling the shear layer thickness, δ. As δ increases, so does the difference between the ideal flux of kinetic energy in the two-layer models and the real flux of kinetic energy in our simulations. A larger difference appears as a larger global

Internal bores: an improved model via a detailed analysis of the energy budget (a) 0.2

301

(b) 0.2 0.1

0 0 –0.2

–0.1 –0.2

–0.4 0

1000

2000

3000

4000

5000

–0.3

0

1

2

3

4

F IGURE 19. Change in the flux of kinetic energy by turbulent mixing in the upper (N), the lower (H), and both ( ) layers as a function of Re (a) and Sc (b) for r = 0.1 and R ≈ 1.9. In the Reynolds-number study, Sc = 1, and in the Schmidt-number study, Re = 3500. The × symbols represent the global energy change due to viscous dissipation.

•

(a)

(b) 0.3 0.4

0.2 0.1

0

0

–0.4

–0.1 –0.8

–0.2 –0.3

–1.2 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8

–0.4

0

0.05

0.10

0.15

0.20

0.25

F IGURE 20. Change in the flux of kinetic energy by turbulent mixing in the upper (N), lower (H), and both ( ) layers as a function of R (a) and r (b) for Re = 3500 and Sc = 1. In the R study, r = 0.1, and in the r study, R ≈ 1.9. The × symbols represent the global energy change due to viscous dissipation.

•

kinetic energy loss, and a larger energy gain in the lower layer, which are consistent with our observation in figure 19(b). As before, the magnitude of energy change due to turbulent mixing is much greater than the change due to viscous dissipation. With regard to the influence of the geometrical parameters, we examine the influence of R in figure 20(a). As R increases, the total amount of kinetic energy loss in a global control volume also increases. The observed change in energy is consistent with the one-layer bore model, which states that the energy loss should increase approximately as R2 . We also observe that the amount of energy gained by the lower layer grows with R until R = 2.4, and then decreases. The trends for r in figure 20(b) are ambiguous. There appears to be a small decrease in the gain of energy across the lower layer as r increases, but the dependence is very weak. As the channel height is decreased for a fixed R, the change in cross-sectional area seen by the upper fluid as it moves across the bore’s front becomes more drastic. Its velocity downstream of the bore must increase. The downstream velocity in the lower layer remains unchanged, so the u-velocity gradient in the z-direction increases

302

Z. Borden, E. Meiburg and G. Constantinescu (b)

0.10 0.05 0

WS

–0.05

KRS

–0.10

(a) 1.6 1.5 1.4 1.3 1.2 1.1 1.0 1.2

1.2

(c)

1.4

1.6

1.8

1.4

1.6

1.8

2.0

2.2

2.4

2.6

2.0

2.2

2.4

2.6

0.8 0.4

1.4

1.6

1.8

2.0

2.2

2.4

2.6

0

R –0.4 –0.8 1.2

R

• •

F IGURE 21. (a) Front velocity as a function of R for both slip ( ) and no-slip (◦) boundaries. (b) Change in energy flux due to viscous forces both globally ( ) and in the lower layer (H). (c) Change in energy flux due to turbulent mixing. Here r = 0.1, Re = 3500 and Sc = 1. Filled symbols represent simulations with slip boundaries and open symbols are for simulations with no-slip boundaries.

near the interface between the two layers, which in turn increases the global energy loss. But increasing r also decreases the front velocity, decreasing the energy loss. The weak dependence of e˙ u and e˙ l on r suggests that these two effects nearly negate each other. In both the R and r parameter studies, the change in kinetic energy flux due to mixing is much greater than the actual energy loss due to viscous dissipation. For geophysically realistic bores, with Re  1, we believe the dominant mechanism for global energy loss and energy gain in the lower layer is turbulent mixing, which changes the balance of kinetic and potential energy in a bore, but is neglected by existing two-layer models. 6.5. Influence of no-slip walls

Earlier, we decided to use free-slip boundaries in our simulations in order to more closely approximate analytical two-layer models. We revisit that decision here to determine if the choice of boundary conditions has a significant impact on a bore’s front velocity or energy flux. Work by H¨artel et al. (2000) found that the bottom boundary condition can have a significant impact on the front velocity and head structure of gravity currents. At a Reynolds number of Re = 4500, the addition of no-slip boundaries decreased the front velocity by 11.1 %. Even at Reynolds numbers as high as Re = 45 000, the front velocity was 5.6 % slower with no-slip boundaries. In contrast to gravity currents, figure 21(a) shows that when we rerun our R parameter study with no-slip boundaries, there is no significant change in the front velocity of internal bores. We therefore expect that the mechanism for energy loss

Internal bores: an improved model via a detailed analysis of the energy budget

303

and transfer within bores should also not be significantly affected by the boundary conditions. The change in energy flux due to turbulent mixing is unaffected by the bottom boundary condition (figure 21c). But in our reference frame moving with the front of the bore, the bottom wall acts as a source of momentum, and viscous effects cause a further increase of energy in the lower layer (figure 21b). For small bores, viscous stresses along the bottom wall cause the total magnitude of the viscous effects to be comparable to the energy change due to mixing. But as can be seen in figure 21(a), the additional energy gain at the bottom boundary does not significantly affect the bore’s height or front velocity. Furthermore, if instead of generating our bores with a dam break, we had generated them with two-layer flow over an obstacle, as has been done in some of the experimental studies of Baines (1984) and Wood & Simpson (1984), the sign of vorticity at the bottom boundary would be opposite, and the no-slip boundary would act as a momentum sink instead of a source. Wood & Simpson (1984) produced laboratory-scale internal bores with both a dam-break and two-layer flow over an obstacle and did not note any differences in the propagation velocities. For these reasons, we believe that only energy changes at the interface between the upper and lower layers affect the bore’s propagation. And in that region, turbulent mixing dominates over viscous effects, regardless of the boundary conditions. 7. Three-dimensional simulations We have used two-dimensional simulations for our parameter and case studies up to this point because they allow us to run simulations at the high grid resolutions we need for direct numerical simulations (DNS) in a reasonable amount of time. But two-dimensional simulations neglect the spanwise instabilities and flow structures that occur in laboratory- and field-scale internal bores. To study these three-dimensional effects and check the validity and generality of our two-dimensional results, we ran some fully three-dimensional simulations. The three-dimensional code we use has been extensively validated (Ooi, Constantinescu & Weber 2007; Gonzalez-Juez et al. 2009) and employs a large-eddy simulation (LES) turbulence model in order to simulate flows at high Reynolds numbers on a coarser grid than would be necessary if were we using DNS. In LES, the Navier–Stokes equations are filtered to eliminate small length scales in the solution. The filtering produces an anisotropic residual stress tensor, which we assume is related to the resolved rate-of-strain tensor through a linear eddy viscosity. With this assumption, our governing equations, (3.1)–(3.3), become

∇ · u = 0,

(7.1)

1 ∂u ∗ + u · ∇u = −∇p + + νSGS ∇ 2 u + ρeg , ∂t Re   ∂ρ 1 ∗ + u · ∇ρ = + κSGS ∇ 2 ρ, ∂t Re Sc 



(7.2) (7.3)

∗ ∗ where νSGS is the dimensionless subgrid-scale viscosity and κSGS is the dimensionless subgrid-scale diffusivity. These two subgrid-scale terms are calculated from the resolved flow field using the dynamic Smagorinsky model (Germano et al. 1991; Lilly 1992). Unfortunately, the code is only second-order accurate in space and time, and as a result, the calculated values of dissipation rates and changes in energy flux will not be as accurate as for the earlier two-dimensional spectral simulations on the

304

Z. Borden, E. Meiburg and G. Constantinescu

(a)

(b) 10 2

3 z 2 1 0

E 0

10 0 10 –2

10

20

x

30

40

50

0

1

2

3 10 –4

y

10 –1

10 0

101

F IGURE 22. (a) Three-dimensional simulation of an internal bore between no-slip walls for Re = 10 000, Sc = 1, R = 2.4 and r = 0.1. The flow field is visualized by the ρ = 0.5 isopycnal and contours of density in the side-plane. (b) One-dimensional turbulence spectra taken on the (x, y)-plane at z = 2.4 showing energy cascades characteristic of fully developed turbulence. Here, E11 refers to the portion of the turbulent kinetic energy resulting from fluctuations in the streamwise velocity. κ1 refers to the wavenumber in the streamwise direction, and κ2 refers to the spanwise wavenumber.

same sized grid. Additional details of the code are provided in Pierce (2001) and Pierce & Moin (2004). The setup for our three-dimensional simulations remains as in figure 4, but now the initial flow field protrudes out of the page in the y-direction. The boundaries are free-slip at the streamwise ends, periodic at the spanwise ends, and no-slip at top and bottom (z = 0 and z = H). We also introduce small, random perturbations into the initial density field to speed the development of spanwise instabilities. We employ the same methods as in the earlier two-dimensional simulations to measure the bore’s front speed and height, but now take an average value in the spanwise direction. The quasi-steady-state density field of an example simulation is shown in figure 22. The evolution of kinetic energy in the LES is given by the following formula adapted from Pope (2000, p. 585):    DEf ∂ 2 r − ρei ui − uj Sij − τij −pδij = Φ − Pr , (7.4) Dt ∂xi Re | {z } |{z} |{z} |{z} (a)

(b)

(c)

(d)

∗ where = is the anisotropic residual stress tensor, and Pr = 2νSGS Sij Sij is the rate of production of residual kinetic energy. We remark that the kinetic energy Ef in (7.4) does not represent the total kinetic energy, but only the kinetic energy of the resolved flow field. For two simulations with different grid sizes, but identical initial conditions, the simulation on the finer grid will contain more resolved kinetic energy after some time t. Furthermore, the difference in the resolved kinetic energy will be at small length scales, and will therefore have a disproportionate impact on the viscous terms in (7.4), which are: (a) the transfer of energy at resolved scales because of viscous stresses; (b) the transfer of energy at resolved scales because of viscous interaction with subgrid-scale modes; (c) the rate of viscous dissipation at the resolved scales; and (d) the rate of energy loss at the resolved scales due to the production of residual kinetic energy. We must choose our grid resolution such that the viscous effects we wish to measure, terms (a,c), are explicitly resolved. The magnitude of (b) should be less than the magnitude of (a), and the magnitude of (d) should be less than the magnitude of

τijr

∗ −2νSGS Sij

Internal bores: an improved model via a detailed analysis of the energy budget (a) 4

(c) 4

3 z 2 1 0

(d) 4

3

3

2

2

z

(b) 4 3 z 2 1 0 –5

2D 3D

305

1 0

5

10

15

20

x

25

30

35

40

0

1

1

0.6 1.0 1.4 1.8

u

F IGURE 23. Density field for R = 2.4, r = 0.1 and Sc = 1 for a two-dimensional (a) and a three-dimensional (b) simulation. The Reynolds numbers are Re = 3500 for (a) and Re = 10 000 for (b). The density field for the three-dimensional simulation is a slice in the (x, z) plane taken at y = 0. (c,d) The mean density and streamwise velocity profiles in the shear layer as functions of z, for both the two-dimensional (2D, dashed line) and three-dimensional (3D, solid line) simulations.

(c), i.e. the LES must effectively be close to a DNS. This constraint limits us to a Reynolds number of Re = 10 000, for a resolution of 3900 × 49 × 256 grid points. Two one-dimensional turbulence spectra, one in the streamwise and one in the spanwise direction, are shown in figure 22(b). Both spectra contain the −5/3 slope characteristic of isotropic turbulence, but only the streamwise spectrum shows that we are resolving the viscous subrange. Our grid in the spanwise direction is too coarse for our simulation to be considered a DNS; however, the spanwise spectrum lacks the pileup of turbulent kinetic energy characteristic of an under-resolved simulation. Our LES model is functioning as it should. We would prefer the grid to be sufficiently resolved in all three dimensions to capture the viscous effects in terms (a,c), but we were unable to refine our grid further. We note, however, that the components of terms (a,c) most responsible for transferring energy between the lower and upper layers are those associated with derivatives in the z-direction. The grid resolution for our simulation is the same in z and x. Since figure 22(b) shows that the grid spacing in x is sufficient to resolve the viscous subrange, we hence expect that the viscous subrange in the wall-normal z-direction is also fully resolved. Therefore, the viscous terms in (7.4) should be representative of the total amount of viscous dissipation and energy transfer. Figure 23, a comparison of our three-dimensional bore to a two-dimensional simulation with the same initial and boundary conditions, shows the density fields of the two simulations at t = 32. The flow structures in the shear layer are more defined in the two-dimensional simulation because the three-dimensional instabilities that cause the breakup of the interfacial vortices are not present. We also notice that the onset of the mixing layer is closer to the front of the current in three dimensions. Figure 23 also shows the average downstream density and streamwise velocity profiles for each of the simulations. The shapes of the profiles are very similar, and all are reasonably well-approximated by a linear stratification. We do, however, observe that the two-dimensional simulations produce slightly thicker shear layers than our three-dimensional simulations. Despite these differences, the measured values of R and uf are similar, with R2D = 2.36, R3D = 2.38, uf ,2D = 1.43 and uf ,3D = 1.44. For comparison, uws = 1.59 and

306 (a)

Z. Borden, E. Meiburg and G. Constantinescu (b) 0.020

0

0.015

–0.02

0.010

–0.04

0.005 –0.06 0 –0.08 –10

(c)

0

10

20

30

40

–10

0

10

20

30

40

0 –0.4 –0.8

Total Lower Upper

–1.2 –1.6 –2.0 –10

0

10

20

30

40

F IGURE 24. Change in the flux of kinetic energy per unit length due to (a) viscous effects, (b) non-hydrostatic effects and (c) turbulent mixing at the interface, plotted as a function of the right-hand boundary of a control volume spanning both layers, or the upper or lower layer separately.

ukrs = 1.49 for the same parameters. The addition of a third dimension increased the front velocity by only 1 %. H¨artel et al. (2000) compared two- and three-dimensional simulations of gravity currents, and also found that large-scale flow properties, such as front velocity and current height, varied by only a few per cent, even though the small-scale flow structures at the interface differed significantly. With (7.4), we attempt to reproduce the results in figures 14–18, where we examined the change in kinetic energy flux in each layer as a function of controlvolume size, for each of our three proposed energy mechanisms. Figure 24 shows that in the three-dimensional simulation, turbulent mixing continues to be the dominant mechanism for global energy loss and lower-layer energy gain. For the present parameter combination, the global effect of turbulent mixing is 33 times greater than the non-hydrostatic effects and 10 times greater than viscous effects. Qualitatively, the effects of non-hydrostatic terms and turbulent mixing are nearly identical for two and three dimensions. The viscous trends, however, differ. Figure 24(a) shows that for the three-dimensional simulation, viscous effects no longer impart energy to the lower layer. The presence of additional spanwise instabilities in three dimensions increases the magnitude of viscous dissipation per unit length, term (c) in (7.4). With term (c), however, the additional spanwise instabilities do not act to diffuse streamwise momentum across the layer interface, only along it. The magnitude of term (a) in (7.4) does not change as much as term (c) on moving to three dimensions. The shift in the relative magnitudes of these two terms means that in three dimensions, the viscous diffusion of momentum into the lower layer is no longer sufficient to overcome the viscous dissipation in that layer, and viscous effects now cause a loss of energy in the lower layer.

Internal bores: an improved model via a detailed analysis of the energy budget

307

Quantitatively, the magnitude of the change in energy flux, per unit length, by turbulent mixing differs modestly between our two- and three-dimensional simulations. Table 1 shows the change in energy flux in each layer, and globally, for three comparable two- and three-dimensional simulations. In all cases, the effect of mixing on the lower layer is reduced in three dimensions. The presence of spanwise instabilities limits the growth of the shear layer’s thickness by suppressing the interfacial vortices, and consequently the amount of mixing is reduced. In other comparisons between two- and three-dimensional simulations, this percentage drop in the effect of mixing lies in a similar range. Therefore, although the magnitude of the effect of mixing is different in two and three dimensions, trends between simulations are unchanged. Based on the results of these three-dimensional simulations, we conclude that our model reduction to two dimensions does not alter our conclusion that turbulent mixing is the dominant mechanism responsible for global energy loss and an energy gain in the lower layer. Furthermore, the trends in e˙ l and e˙ u as we vary parameters are preserved when moving from two- to three-dimensional simulations. The magnitudes of the terms, however, are reduced. 8. Improved internal bore model Turbulent mixing at the interface between the upper and lower layers causes a redistribution of kinetic energy that is not accounted for by the control-volume models of internal bores. We incorporate the effects of this mixing into the LC relation by comparing the flux of kinetic energy downstream of the bore for the ‘ideal’ and ‘real’ velocity profiles in figure 9. We assume that near the interface, the horizontal velocity transitions linearly between U1 and U2 over some distance δ. Our assumption of a linear velocity profile is common in studies of stratified shear layers (Turner 1986) and represents a good approximation to the real velocity profile, seen in figure 23. The flux of kinetic energy in the ‘ideal’ and ‘real’ cases will be identical for all z except (hf − δ/2) < z < (hf + δ/2). If we define z∗ = 0 at the bottom of the shear layer, the ideal and real kinetic energy fluxes in the shear layer region of the lower layer downstream of the bore are Z 1 δ/2 3 ∗ U1 dz , (8.1) e˙ ideal = 2 0  3 Z 1 δ/2 z∗ e˙ real = U1 + (U2 − U1 ) dz∗ . (8.2) 2 0 δ

The apparent energy change in the lower layer is the difference between these two terms. Integrating and taking the difference, we find that  δ  24U12 (U2 − U1 ) + 8U1 (U2 − U1 )2 + (U2 − U1 )3 . (8.3) 128 We can now use (6.8) to relate the thickness of the shear layer to the change in velocity across the layers. We remark that Turner (1986) noticed large discrepancies between the thickness values predicted by (6.8) and experimental observations. He argued that viscous energy losses in the shear layer cause (6.8) to overpredict the thickness of the shear layer and that his expression should be treated as a scaling argument rather than an exact expression. If we apply his result to (8.3), and write e˙ l =

e˙ l e˙ u e˙ t e˙ l e˙ u e˙ t e˙ l e˙ u e˙ t

3.23 × 10−3 −5.11 × 10−3 −1.89 × 10−3 −2.09 × 10−2 −6.34 × 10−2 −8.43 × 10−2 −1.83 × 10−2 −4.70 × 10−2 −6.53 × 10−2

3.00 × 10−4 −2.73 × 10−2 −2.70 × 10−2 1.35 × 10−1 −6.23 × 10−1 −4.88 × 10−1 6.85 × 10−2 −5.75 × 10−1 −5.07 × 10−1

Three dimensional Viscous Mixing 7.37 × 10−3 −3.33 × 10−3 −1.07 × 10−2 −1.64 × 10−3 −5.95 × 10−2 −6.11 × 10−2 −1.64 × 10−3 −5.95 × 10−2 −6.11 × 10−2

3.64 × 10−2 −1.58 × 10−2 −5.25 × 10−2 5.52 × 10−1 −1.06 × 100 −5.08 × 10−1 5.52 × 10−1 −1.06 × 100 −5.08 × 10−1

Two-dimensional Viscous Mixing −56.19 53.53 −82.83 1177.34 6.59 38.03 1019.73 −21.05 6.89

−99.18 73.48 −48.48 −75.53 −41.22 −3.99 −87.57 −45.75 −0.36

% change from 2D Viscous (%) Mixing (%)

TABLE 1. Change in energy flux due to viscous effects and mixing for two- and three-dimensional simulations. In trial 1, R = 1.37, r = 0.1, Re = 10 000 and Sc = 1. In trial 2, R = 2.40, r = 0.1, Re = 5000 and Sc = 1. In trial 3, R = 2.40, r = 0.1, Re = 10 000 and Sc = 1.

3

2

1

Trial

308 Z. Borden, E. Meiburg and G. Constantinescu

Internal bores: an improved model via a detailed analysis of the energy budget

309

the resulting expression as a function of the bore’s front velocity, we arrive at the following expression for the energy change in the lower layer: e˙ l = where

β(R − 1)u5f R2 (1 − Rr)

,

(8.4)

 (R − 1)2  2 2 24 (1 − Rr) +8(R − 1)(1 − Rr) + (R − 1) , (8.5) R3 (1 − Rr)4 and f is Turner’s scaling coefficient, which includes the functional dependence of δ on the Reynolds and Schmidt numbers. Koop & Browand (1979), Turner (1986), and our own observations of energy gain in the lower layer suggest that δ does not depend strongly on the Reynolds number as long as it is sufficiently high. We do, however, observe a strong dependence of δ on mass diffusion. As a first approximation, we set f = CPe−1/2 , where C is some constant, and Pe = Re Sc is the P´eclet number. This scaling is appropriate for a laminar concentration boundary layer. To derive the proper expression for the front velocity, we combine (8.4) with the LC relation as stated in (2.20)
1 R2 r − 3Rr + R + 1 u2f = − 4β i1/2 1 h 2 2 + (R r − 3Rr + R + 1) −8βR2 (1 − Rr)(Rr + r − 2) . (8.6) 4β β =f

First, we compare (8.6), our new expression for uf , with the results of our parameter studies, and plot the results in figure 25. Our scaling analysis, when substituted into the LC relation, accurately predicts the front velocity of our simulated bores over a wide range of parameter values. The best fit is obtained with C = 0.5. Second, we examine how well our scaling analysis predicts the energy gain in the lower layer, keeping in mind that in § 7 we found that the absolute, but not the relative, magnitudes of the rate of change of energy can be different between two- and three-dimensional simulations. Figure 26 shows that (8.4) makes reasonably accurate predictions with two exceptions. We do not capture the correct energy gain at low Reynolds numbers when viscous effects, which are ignored by our model, become important. We also underpredict the gain of energy in the lower layer, and overpredict the front velocity, for small values of R. Benjamin & Lighthill (1954) proposed that undular bores, which arise for small values of R, are capable of radiating energy through the production of internal waves at their fronts. We do not capture this energy transport mechanism, and hence we cannot accurately model undular bores. But for these cases, the present model is no worse than either the KRS or WS models, just not much better. Another shortcoming is that our model does not predict the non-monotonic scaling of e˙ l on R. Our model predicts that e˙ l should always increase with increasing R, instead of an increase followed by a decrease. In order to model large bores accurately enough to reproduce this trend, we would need to include the presence of the shear layer in our statements of conservation of mass and momentum because in large bores, the shear layer is very thick compared to ha . Klemp et al. (1997) derived such a model by assuming an ‘integrated velocity deficit’ in each layer, which they termed Q1 and Q2 . They included these velocity deficits in the equations for conservation of mass and momentum, and assumed that the velocity deficit becomes sufficiently small

310

Z. Borden, E. Meiburg and G. Constantinescu 1.4 1.3 1.2 1.1

1.6 0

1000

2000

3000

4000

5000

1.4 1.2 1.0

1.40

1.5

2.0

2.5

1.35 1.30 1.25 1.20

1.6 0

1000

2000

3000

4000

5000 1.4 1.2

1.40

1.0

1.35

0

1.30

0.05

0.10

0.15

0.20

0.25

1.25 1.20

0

1

2

3

4

•

F IGURE 25. Comparison of the the front velocities in our parameter studies (, , H) with WS model (dotted lines), the KRS model (dashed lines), and our new model given by (8.6) (solid lines). The parameter values in all simulations, unless otherwise stated, are R = 1.9, r = 0.1, Sc = 1 and Re = 3500.

behind the bore such that high-order terms can be neglected. But their solution was still underdetermined – a further relation is needed to solve for the ratio Q1 /Q2 . This model does offer an advantage over ours in that it includes the effects of mixing on mass and momentum conservation. But it is unclear whether or not their approach would be any more accurate for larger bores because of their assumption of a small downstream velocity deficit. Third, we examine the total loss of energy across a bore. To calculate e˙ t , we take our front velocity from (8.6), plug it into (2.19) along with R and r, and add the resulting value of e˙ u to the value of e˙ l computed by (8.4). Figure 27 shows that although our prediction of the total energy loss across a bore is not as good as our prediction of the lower-layer energy gain, it is still a reasonable estimation of the energy loss across a bore due to mixing. Finally, we note that although this model provides a more accurate description of the front velocity than both the WS and KRS models, as well as yielding an estimate of the amount of energy lost to mixing, it does not illuminate the reason why a two-layer, inviscid control-volume method would predict an energy loss in the first place. When first deriving his jump condition for gravity currents, Benjamin (1968) noticed that only currents whose height equals half the channel height can be energy conserving. A shallower current, as would be produced by a partial-depth lock release, has to be accompanied by a loss of energy. Similarly, in the two-layer bore models, the only

Internal bores: an improved model via a detailed analysis of the energy budget

311

0.25 0.20 0.15 0.10 0.05 0

1000

2000

3000

4000

5000

0.20 0.15

0.5 0.4 0.3 0.2 0.1 0 1.0

1.5

2.0

2.5

3.0

0.10 0.05 0.20 0

1000

2000

3000

4000

5000

0.15 0.10

0.20

0.05

0.15

0

0.10

0.05

0.10

0.15

0.20

0.25

0.05 0

1

2

3

4

F IGURE 26. Comparison of the energy change in the lower layer given by our parameter studies (, , H) to the value predicted by our model for e˙ l (solid lines). The parameter values in all simulations, unless otherwise stated, are R = 1.9, r = 0.1, Sc = 1 and Re = 3500.

•

way a bore can have no loss of energy is if hf = H/2. If the height is any less, the two-layer models predict a global loss of energy. The underlying reason must be a more fundamental problem with the control-volume models, related to the fact that the calculations of energy flux across an internal bore use the square of a depth-averaged velocity, not the depth average of the velocity squared. 9. Conclusions The present investigation employs two- and three-dimensional highly resolved numerical simulations to evaluate two-layer hydrostatic models for internal bores. The main objective was to search for a mechanism whereby the lower, expanding layer of an internal bore could be gaining energy across the bore. An energy gain in the lower layer would explain why the KRS model, despite being a much better fit for the experimental data than the WS model, still overpredicts the front velocity for certain parameter ranges. We studied a representative bore in detail and found that turbulent mixing at the interface between the two layers, which is neglected in the two-layer control-volume models, converts some of the kinetic energy into potential energy, and makes it appear as if the lower layer has gained energy. We then carried out a parameter study and proposed a new model for internal bores, which incorporates the difference in kinetic energy flux between an ideal velocity field and a velocity field with a linear shear layer as a gain of energy by the lower layer. We obtain

312

Z. Borden, E. Meiburg and G. Constantinescu 0 –0.1 –0.2 –0.3 –0.4

0

1000

2000

3000

4000

5000

0

0 –0.2 –0.4 –0.6 –0.8 –1.0 1.0

–0.1

1.5

2.0

2.5

3.0

–0.2 –0.3 –0.4

0 0

1000

2000

3000

4000

5000

–0.05 –0.10

0

–0.15

–0.05

–0.20

–0.10

–0.25

0

0.05

0.10

0.15

0.20

0.25

–0.15 –0.20 –0.25

0

1

2

3

4

•

F IGURE 27. Comparison of the global energy change given by our parameter studies (, , H) to the value predicted by our model for e˙ t and the LC relation (solid lines). The parameter values in all simulations, unless otherwise stated, are R = 1.9, r = 0.1, Sc = 1 and Re = 3500.

a closed-form solution from the LC relation that yields improved agreement with experimental data, compared to the WS and KRS models. An assumption made throughout this investigation was that our internal bores propagated along an interface between two miscible fluids. While it is certainly the case that most naturally occurring Boussinesq internal bores propagate along a miscible interface, it is still interesting to consider how having an immiscible interface would affect the bore propagation. We suspect that in flows with a high capillary number, where viscous forces dominate over surface tension, mixing and entrainment will still be prominent at the interface. These bores should propagate similarly to miscible ones. For bores with low capillary number, however, mixing will definitely be suppressed by surface tension forces and our model would no longer be applicable. We remark that the results of this investigation are of interest not only for internal bores. Analytical models of gravity currents also predict a loss of energy across the front of a gravity current, as long as the current occupies less than half of the channel depth. We would not be surprised if turbulent mixing and entrainment is the responsible mechanism in gravity currents as well. On a final note, we comment that Li & Cummins (1998) used experimental data from Baines (1984) to argue against a gain of energy in the lower layer. They found that most of Baines’ experimental front velocities fell in between the limits given by the WS and KRS models. They therefore reason that both layers must have a net loss

Internal bores: an improved model via a detailed analysis of the energy budget

313

of energy. Baines’ experimental bores, however, propagated along an interface between kerosene and water, two immiscible fluids with a density ratio ρ2 /ρ1 = 0.82. For this density ratio, non-Boussinesq effects may become important. Hence, it will be of interest to investigate the extent to which the present findings apply to non-Boussinesq bores. Acknowledgements This research has been supported by NSF grants CBET-0854338, CBET-1067847 and OCE-1061300. Our simulations were carried out at the Community Surface Dynamics Modelling System (CSDMS) high-performance computing facility at the University of Colorado in Boulder. We thank CSDMS director Professor J. Syvitski and the technical staff at CSDMS for their support. Finally, we thank Dr J. Rottman at the University of California, San Diego for bringing this problem to our attention and for his help throughout this work. REFERENCES

BAINES, P. G. 1984 A unified description of two-layer flow over topography. J. Fluid Mech. 146, 127–167. BAINES, P. G. 1995 Topographic Effects in Stratified Fluids. Cambridge University Press. B ENJAMIN, T. B. 1968 Gravity currents and related phenomena. J. Fluid Mech. 31, 209–248. B ENJAMIN, T. B. & L IGHTHILL, M. J. 1954 On cnoidal waves and bores. Proc. R. Soc. Lond. A 224 (1159), 448–460. B ONOMETTI, T. & BALACHANDAR, S. 2008 Effect of Schmidt number on the structure and propagation of density currents. Theor. Comput. Fluid Dyn. 22, 341–361. C HU, V. H. & BADDOUR, R. E. 1977 Surges, waves and mixing in two-layer density stratified flow. In Proc. 17th Congress of the International Association of Hydraulic Research, Baden-Baden, Germany, pp. 303–310. IAHR. C LARKE, R. H., S MITH, R. K. & R EID, D. G. 1981 The morning glory of the Gulf of Carpentaria: an atmospheric undular bore. Mon. Weath. Rev. 109, 1726–1750. C UMMINS, P. F. 1995 Numerical simulations of upstream bores and solitons in a two-layer flow past an obstacle. J. Phys. Oceanogr. 25, 1504–1515. F LETCHER, C. A. J. 1991 Computational Techniques for Fluid Dynamics, vol. 2, 2nd edn. Springer. G ERMANO, M., P IOMELLI, U., M OIN, P. & C ABOT, W. H. 1991 A dynamic subgrid-scale eddy viscosity model. Phys. Fluids 3 (7), 1760–1765. G ONZALEZ -J UEZ, E. & M EIBURG, E. 2009 Shallow-water analysis of gravity-current flows past isolated obstacles. J. Fluid Mech. 635, 415–438. G ONZALEZ -J UEZ, E., M EIBURG, E. & C ONSTANTINESCU, G. 2009 Gravity currents impinging on bottom-mounted square cylinders: flow fields and associated forces. J. Fluid Mech. 631, 65–102. G ONZALEZ -J UEZ, E., M EIBURG, E., T OKYAY, T. & C ONSTANTINESCU, G. 2010 Gravity current flow past a circular cylinder: forces, wall shear stresses and implications for scour. J. Fluid Mech. 649, 69–102. H AASE, S. P. & S MITH, R. K. 1984 Morning glory wave clouds in Oklahoma: a case study. Mon. Weather Rev. 112, 2078–2089. ¨ H ARTEL , C., M EIBURG, E. & N ECKER, F. 2000 Analysis and direct numerical simulation of the flow at a gravity-current head. Part 1. Flow topology and front speed for slip and no-slip boundaries. J. Fluid Mech. 418, 189–212. H ORN, D. A., I MBERGER, J. & I VEY, G. N. 2001 The degeneration of large-scale interfacial gravity waves in lakes. J. Fluid Mech. 434, 181–207. H OSEGOOD, P. & VAN H AREN, H. 2004 Near-bed solibores over the continental slope in the Faeroe–Shetland channel. Deep-Sea Res. II 51 (25–26), 2943–2971.

314

Z. Borden, E. Meiburg and G. Constantinescu

H UPPERT, H. E. & S IMPSON, J. E. 1980 The slumping of gravity currents. J. Fluid Mech. 99, 785–799. K LEMP, J. B., ROTUNNO, R. & S KAMAROCK, W. C. 1994 On the dynamics of gravity currents in a channel. J. Fluid Mech. 269, 169–198. K LEMP, J. B., ROTUNNO, R. & S KAMAROCK, W. C. 1997 On the propagation of internal bores. J. Fluid Mech. 331, 81–106. KOOP, C. G. & B ROWAND, F. K. 1979 Instability and turbulence in a stratified fluid with shear. J. Fluid Mech. 93 (1), 135–179. L ANDAU, L. D. & L IFSHITZ, E. M. 1959 Fluid Mechanics. Pergamon. L EICHTER, J. J., W ING, S. R., M ILLER, S. L. & D ENNY, M. W. 1996 Pulsed delivery of subthermocline water to Conch Reef (Florida Keys) by internal tidal bores. Limnol. Oceanogr. 41 (7), 1490–1501. L I, M. & C UMMINS, P. F. 1998 A note on hydraulic theory of internal bores. Dyn. Atmos. Oceans 28, 1–7. L ILLY, D. K. 1992 A proposed modification of the Germano-subgrid-scale closure method. Phys. Fluids A 4 (3), 633–635. M OROZOV, E. G., T RULSEN, K., V ELARDE, M. G. & V LASENKO, V. I. 2002 Internal tides in the straight of gibraltar. J. Phys. Ocenaogr. 32, 3193–3206. ¨ N ECKER, F., H ARTEL , C., K LEISER, L. & M EIBURG, E. 2005 Mixing and dissipation in particle-driven gravity currents. J. Fluid Mech. 545, 339–372. O OI, S. K., C ONSTANTINESCU, S. G. & W EBER, L. 2007 Numerical simulations of lock-exchange compositional gravity currents. J. Fluid Mech. 635, 361–388. PANTON, R. L. 2005 Incompressible Flow, 3rd edn. John Wiley. P IERCE, C. D. 2001 Progress-variable approach for large eddy simulation of turbulent combustion. PhD thesis, Stanford University. P IERCE, C. D. & M OIN, P. 2004 Progress-variable approach for large-eddy simulation of non-premixed turbulent combustion. J. Fluid Mech. 504, 73–97. P OPE, S. B. 2000 Turbulent Flows. Cambridge University Press. S IMPSON, J. E. 1997 Gravity Currents, 2nd edn. Cambridge University Press. T URNER, J. S. 1973 Buoyancy Effects in Fluids. Cambridge University Press. T URNER, J. S. 1986 Turbulent entrainment: the development of the entrainment assumption, and its application to geophysical flows. J. Fluid Mech. 173, 431–471. U NGARISH, M. 2008 Energy balances and front speed conditions of two-layer models for gravity currents produced by lock release. Acta Mech. 201, 63–81. U NGARISH, M. 2009 Energy balances for gravity currents with a jump at the interface produced by lock release. Acta Mechanica 211, 1–21. U NGARISH, M. & H UPPERT, H. E. 2006 Energy balances for propagating gravity currents: homogeneous and stratified ambients. J. Fluid Mech. 565, 363–380. U NGARISH, M. & H UPPERT, H. E. 2008 Energy balances for axisymmetric gravity currents in homogeneous and linearly stratified ambients. J. Fluid Mech. 616, 303–326. WAKIMOTO, R. M. & K INGSMILL, D. E. 1995 Structure of an atmospheric undular bore generated from colliding boundaries during CaPE. Mon. Weather Rev. 123, 1374–1393. W INANT, C. D. & B ROWAND, F. K. 1974 Vortex pairing: the mechanism of turbulent mixing-layer growth at moderate Reynolds number. J. Fluid Mech. 63 (2), 237–255. W OOD, I. R. & S IMPSON, J. E. 1984 Jumps in layered miscible fluids. J. Fluid Mech. 140, 215–231. Y IH, C. S. & G UHA, C. R. 1955 Hydraulic jump in a fluid system of two layers. Tellus 7, 358–366. Z HOU, J., D UPUY, B., B ERTOZZI, A. L. & H OSOI, A. E. 2005 Theory for shock dynamics in particle–laden thin films. Phys. Rev. Lett. 94, 117803.