2.1. Problem description

We consider an internal hydraulic jump propagating into a quiescent ambient fluid, generated, for example, by a dam-break release of an initially discontinuous interface (see figure 1). Through the jump, the density ${\it\rho}(z)$ and isopycnals (streamlines) are displaced a distance ${\it\eta}(z)$ relative to the upstream ambient, where the density profile is ${\it\rho}_{a}(z-{\it\eta})$ . The velocity is $U(z)$ in a frame of reference moving with speed $C_{b}$ , which is assumed constant (the velocity is positive to the right). The isopycnal (streamline) displacement relative to the ambient is  ${\it\eta}(z)$ .

Figure 1.

Schematic for an internal bore in continuous stratification. Velocity and density profiles are $U(z)$ and ${\it\rho}(z)$ , where $z$ is the vertical coordinate through the bore region. The isopycnal displacement, ${\it\eta}(z)$ , is measured relative to the upstream ambient, where the density profile is ${\it\rho}_{a}(z-{\it\eta})$ . The bore speed is $C_{b}$ .

2.2. Internal bores in two-layer stratification

Several models have been proposed for the two-layer limit that predict the bore speed, $C_{b}$ , as a function of the interface location in the ambient, $h_{a}$ (measured from the bottom), and the bore height, $h_{b}$ . To obtain jump conditions, mass and horizontal momentum conservation laws are enforced in each layer separately. However, an additional assumption is required to obtain an expression for the pressure jump across the bore. Several closures have been proposed, based on assumptions about the distribution of pressure (Yih & Guha Reference Yih and Guha1955), or more often on the energy dissipation across the jump (see Baines Reference Baines1995). The theory by Chu & Baddour (Reference Chu and Baddour1977) and Wood & Simpson (Reference Wood and Simpson1984) (often termed CBWS) assumes that energy is conserved in the contracting layer and all dissipation occurs in the expanding layer. Klemp, Rotunno & Skamarock (Reference Klemp, Rotunno and Skamarock1997) (termed KRS) developed a closure that conserves energy in the expanding layer. The KRS closure reproduces the Benjamin (Reference Benjamin1968) gravity current front condition in the limit $h_{a}\rightarrow 0$ , while that of CBWS does not. Li & Cummins (Reference Li and Cummins1998) (termed LC) detailed a more general theory in which the energy loss could be distributed arbitrarily between layers, yielding a family of solutions for the bore speed.

Borden, Meiburg & Constantinescu (Reference Borden, Meiburg and Constantinescu2012), in a numerical investigation of the energy budgets for two-layer bores, found that energy could be transferred from the upper to the lower layer by turbulent mixing. They derived a closure for the bore speed based on a semi-empirical relationship for the shear-layer thickness associated with turbulent mixing and the associated energy dissipation, termed the Borden–Meiburg circulation (BMC) model. In subsequent work (Borden & Meiburg Reference Borden and Meiburg2013b ) the same authors introduced a model to parametrize the turbulent mixing across the bore, arguing that conservation of vorticity (rather than energy) should be enforced through the bore to account for the fact that vertical momentum is not otherwise considered in the conservation laws. This resulted in two new closures, one in the limit where the sharp two-layer stratification is preserved through the bore, termed the vortex sheet (VS) model, and an additional closure that accounts for interface thickening due to turbulent mixing, termed the diffuse vortex sheet (DVS) model.

Despite the ubiquity of two-layer theories, the problem continues to be open. This is in part because none of these theories has yet produced ideal agreement with observations and numerical results over the full range of parameter values. One major issue is that they all assume steady flow in a frame moving with the bore. This is an inherent simplification, because internal hydraulic jumps are in fact often dispersive. Recent work by Esler & Pearce (Reference Esler and Pearce2011) investigated unsteady undular bores using an extension of Whitham (Reference Whitham1974) modulation theory applied to a fully nonlinear, weakly dispersive internal wave model.

The inclusion of downstream interface thickness in the BMC and DVS models and their apparent improvement over sharp two-layer models suggests that continuous stratification is an important consideration. The approach we describe incorporates stratification through first principles by directly considering the ambient density structure in the conservation laws. As a result, we are able to predict the velocity and density structure through the bore a priori, resulting in improved agreement with simulation results without the need for semi-empirical mixing parametrizations. Moreover, our results suggest that it is in fact the ambient waveguide, rather than turbulent mixing, that structures internal bore propagation. We return to a discussion of two-layer theories in § 3 after first discussing our continuous model.

2.3. Nonlinear internal waves and conjugate states in continuous stratification

To develop our theory, we view an internal bore as a limiting form of a nonlinear solitary wave, which can be described by the DJL equation (Dubreil-Jacotin Reference Dubreil-Jacotin1934; Long Reference Long1953; Stastna & Lamb Reference Stastna and Lamb2002). In a frame of reference moving with the solitary-wave speed $c$ , the streamline displacement, ${\it\eta}(x,z)$ , measured relative to the ambient with density profile ${\it\rho}_{a}(z_{a})$ , where $z_{a}=z-{\it\eta}(x,z)$ , is given by

subject to boundary conditions ${\it\eta}=0$ at $z=0,H$ and ${\it\eta}\rightarrow 0$ as $x\rightarrow \pm \infty$ . Here $N$ is the Brunt–Väisälä frequency, $N^{2}=-(g/{\it\rho}_{o})\,\text{d}{\it\rho}_{a}/\text{d}z_{a}$ , and ${\it\rho}_{o}$ is a constant reference density. While (2.1) neglects ambient shear, it can be incorporated (see Stastna & Lamb Reference Stastna and Lamb2002).

Lamb (Reference Lamb2002) showed that the wave amplitude increases with $c$ up to one of two limiting outcomes (depending on ambient stratification and shear). (i) Overturning occurs when the local velocity matches the wave speed ( $u=c$ in the laboratory frame). A necessary condition for overturning is that $N$ be non-zero at the boundary, i.e.  $z=0\ (z=1)$ for waves of elevation (depression) (Lamb Reference Lamb2002). (ii) DJL solutions become broad-crested and approach an infinite-wavelength solution that, in the half-plane, represents an energy-conserving bore (Brown & Christie Reference Brown and Christie1998; Lamb & Wan Reference Lamb and Wan1998). In this limit, the flat isopycnals through the wave are said to be conjugate to those in the ambient, and we refer to this solution as the conjugate state. Lamb & Wan (Reference Lamb and Wan1998) found conjugate-state solutions for arbitrary stratification and shear. In the two-layer Boussinesq limit, Lamb (Reference Lamb2000) showed that the conjugate-state speed is $C_{cs}=0.5(g^{\prime }H)^{1/2}$ , and the upstream pycnocline depth is found at mid-depth $h_{cs}=H/2$ , both independent of $h_{a}$ , recovering the energy-conserving CBWS and KRS two-layer solutions.

3.1. Conservation equations

The conservation equations for an internal jump in continuous stratification follow the development in Lamb & Wan (Reference Lamb and Wan1998), Lamb & Wilkie (Reference Lamb and Wilkie2004) and White & Helfrich (Reference White and Helfrich2008). Referring to figure 1, we assume that there exists a region behind the bore in which the flow is uniform in $x$ , and the pressure hydrostatic. Considering nominally inviscid motion slightly modified by dissipation, Bernoulli’s equation can be written in a frame moving with $C_{b}$ , along a streamline between the ambient (where variables have subscript $a$ ) and the region behind the bore (without subscripts) as

Here ${\it\Delta}(z)$ represents a Bernoulli head loss with arbitrary vertical structure. Here and subsequently, we make the Boussinesq approximation by assuming that density variation is weak relative to a fixed reference value, ${\it\rho}_{o}$ , and is therefore only important in the gravity term. Diffusion of density is also neglected, implying that density isopycnals coincide with streamlines, i.e.  ${\it\rho}={\it\rho}({\it\psi})$ . By continuity, the velocity in the moving frame is related to the bore speed by $U(z)=C_{b}({\it\eta}_{z}-1)$ (positive velocity is in the direction of bore propagation). Hydrostatic pressure gives $\text{d}p/\text{d}z=-{\it\rho}_{a}(z-{\it\eta})g$ and $\text{d}p_{a}/\text{d}z_{a}=-{\it\rho}_{a}(z_{a})g$ through the bore and ambient regions, respectively. After (i) taking the vertical derivative, $\text{d}/\text{d}z$ , of (3.1), (ii) substituting the continuity and hydrostatic pressure relations, and (iii) noting that a vertical derivative of any quantity is related to its derivative in the ambient by $\text{d}/\text{d}z=(1-{\it\eta}_{z})\,\text{d}/\text{d}z_{a}$ , an equation is obtained for the streamline displacement,

which is subject to the boundary conditions ${\it\eta}(0)={\it\eta}(H)=0$ . This constitutes an eigenvalue problem for the bore speed, $C_{b}$ , and the streamline displacement, ${\it\eta}(z)$ . In the linear limit ${\it\eta}\rightarrow 0$ and with ${\it\Delta}=0$ , (3.2) reduces to the linear equation governing normal modes in a continuously stratified Boussinesq fluid.

As with two-layer theories, conservation of momentum must also be enforced between the bore and ambient regions, resulting in an additional constraint,

Using (3.1) to write $p_{a}(z_{a})$ in terms of $p(z)$ and again using $\text{d}/\text{d}z=(1-{\it\eta}_{z})\,\text{d}/\text{d}z_{a}$ , the hydrostatic relation and integration by parts yields

In the special case ${\it\Delta}=0$ , (3.2) and (3.4) are equivalent to the conjugate state of Lamb & Wan (Reference Lamb and Wan1998). The addition of the arbitrary ${\it\Delta}(z)$ term allows for a range of solutions, but requires an appropriate closure. Complete derivations of the non-Boussinesq versions of both (3.2) and (3.4) are given in appendix A.

3.2. Bernoulli head loss

We assume that the head loss can be written as the product of a constant and a vertical shape function, ${\it\Delta}(z)={\it\Delta}_{o}\,f(z)$ . Given an assumption for $f(z)$ , ${\it\Delta}_{o}$ (along with $C_{b}$ ) is part of the solution obtained from (3.2) and (3.4). The shape function distributes the head loss across isopycnals, and is the continuous analogue of the two-layer closures that distribute dissipation arbitrarily between layers.

To illustrate, consider the following distribution, which is linear in the non-dimensional density, or buoyancy, $b(z)=({\it\rho}-{\it\rho}(H))/({\it\rho}(0)-{\it\rho}(H))$ ,

where ${\it\Delta}_{o}={\it\Delta}(0)+{\it\Delta}(H)$ and ${\it\epsilon}=({\it\Delta}(0)-{\it\Delta}(H))/{\it\Delta}_{o}$ . Here ${\it\Delta}_{o}/2$ is the arithmetic mean head loss (and the loss on the mid-density isopycnal, ${\it\Delta}_{o}/2={\it\Delta}_{b=0.5}$ ), while ${\it\epsilon}$ measures the difference relative to the mean. The one-parameter family in ${\it\epsilon}$ distributes dissipative losses non-uniformly across isopycnals, and recovers the LC general theory in the two-layer limit, as illustrated by the following cases.

1. (a) ${\it\epsilon}=1$ . Here ${\it\Delta}(z)={\it\Delta}_{o}b(z)$ , confining the loss near the bottom expanding region. In the two-layer limit, $b(z)$ is a step function and (3.5) is equivalent to the CBWS closure.

2. (b) ${\it\epsilon}=-1$ . Here ${\it\Delta}={\it\Delta}_{o}(1-b(z))$ , which distributes the head loss through the upper contracting region. In the two-layer limit, (3.5) is equivalent to the KRS solution.

3. (c) ${\it\epsilon}=0$ . Here ${\it\Delta}={\it\Delta}_{o}/2=\text{const.}$ , which distributes the head loss uniformly. In the two-layer limit, this is equivalent to the Borden & Meiburg (Reference Borden and Meiburg2013b ) vorticity conservation VS model (see § 3.3).

4. (d) ${\it\epsilon}\rightarrow -\infty$ . In this case, ${\it\Delta}(z)={\it\Delta}_{o}{\it\epsilon}(b(z)-(1/2))$ and ${\it\Delta}(0)=-{\it\Delta}(H)$ , implying that energy is gained in the lower layer at the expense of the upper layer. Although ${\it\Delta}_{o}\rightarrow 0$ , the product ${\it\epsilon}{\it\Delta}_{o}$ is finite, and it will be shown later that this solution results, perhaps counter-intuitively, in a larger rate of total dissipation than the other special cases (a)–(c).

In the two-layer limit, the general LC theory can be written in terms of ${\it\epsilon}$ as

where $R=h_{b}/h_{a}$ is the ratio of the bore height to the ambient interface height, $r=h_{a}/H$ and $g^{\prime }=g({\it\rho}(0)-{\it\rho}(H))/{\it\rho}_{o}$ . It can be verified that (3.6) recovers the CBWS and KRS solutions for ${\it\epsilon}=1$ and ${\it\epsilon}=-1$ , respectively. For continuous stratification, the structure of ${\it\Delta}(z)$ may be more complex. In § 5.5 we compare (3.5) with profiles of dissipative losses from Navier–Stokes simulations.

With the head loss given by (3.5), the eigenvalue problem governing the streamline displacement through the bore (3.2) reduces to

3.3. Vorticity production

If (3.2) is multiplied by $C_{b}$ , one then obtains a general statement about vorticity production downstream of the bore,

Downstream vorticity arises from baroclinicity (first term on right) and gradients of dissipation across streamlines (lines of constant $C_{b}z_{a}$ ) (second term on right). In the absence of baroclinic effects, this result has a well-known counterpart in rotating single-layer shallow-water flows, where changes in potential vorticity $q$ across a steady jump are related to gradients of the Bernoulli function across streamlines through $[q]=\text{d}[B]/\text{d}{\it\psi}$ , where $B$ is the Bernoulli function and ${\it\psi}$ is the streamfunction (see Pratt & Whitehead Reference Pratt and Whitehead2008, § 3.5). When ${\it\epsilon}=0$ in (3.5), the dissipation is the same along all streamlines, and downstream vorticity arises solely from baroclinic effects. Thus the two-layer VS model of Borden & Meiburg (Reference Borden and Meiburg2013b ), which is based on vorticity production by baroclinic effects only, is recovered from (3.6) with ${\it\epsilon}=0$ . (This result is also noted more recently by Borden & Meiburg (Reference Borden and Meiburg2013a ).)

3.4. Energy dissipation

The total rate of energy dissipation for a steady Boussinesq stratified bore is given by the difference in energy flux between sections b and a, respectively, in figure 1,

Through manipulations similar to those described in § 3.1, (3.9) can be simplified to

where $D>0$ implies a net energy loss, as required for a physically realistic bore solution.

In the two-layer limit, this expression reduces to $D=e_{l}+e_{u}$ , the sum of the dissipation in the lower and upper layers, respectively, as given by the formulas (9) and (11) in Li & Cummins (Reference Li and Cummins1998). In terms of ${\it\Delta}_{o}$ and ${\it\epsilon}$ ,

These expressions illustrate that the case with ${\it\epsilon}\rightarrow -\infty$ and ${\it\Delta}_{o}{\it\epsilon}$ finite, i.e.  $D=-{\it\Delta}_{o}{\it\epsilon}C_{b}H(1-2r)/2$ , produces the maximum dissipation compared with other ${\it\epsilon}$ values. Even though the head losses are equal and opposite, ${\it\Delta}(0)=-{\it\Delta}(H)={\it\Delta}_{o}{\it\epsilon}/2$ , the loss in the contracting layer is carried by a faster current, $C_{b}(1-r)/(1-Rr)$ , so that it more than offsets the energy gain in the lower layer with slower speed, $C_{b}/R$ , such that $|e_{u}|>|e_{l}|$ and $D>0$ .

3.5. Trapped cores

It is possible for (3.2) to develop solutions for which $z_{a}=z-{\it\eta}(z)<0$ , which implies, unphysically, that isopycnals originate from $z_{a}<0$ (Lamb & Wan Reference Lamb and Wan1998; Stastna & Lamb Reference Stastna and Lamb2002; Helfrich & White Reference Helfrich and White2010). The onset of this condition occurs when ${\it\eta}_{z}=1$ , or $U(z)=0$ in the wave frame. Lamb & Wilkie (Reference Lamb and Wilkie2004) found solutions to the conjugate DJL equation beyond this limit by assuming a uniform trapped core of constant density ${\it\rho}={\it\rho}_{a}(0)$ (for bottom-propagating waves), and Helfrich & White (Reference Helfrich and White2010) used a similar framework to find solutions to the full DJL equation (2.1) with a trapped core. In these methods, the flow must be matched along the core boundary with the DJL equation outside. A simpler approach, suggested by Helfrich & White (Reference Helfrich and White2010) and King, Carr & Dritschel (Reference King, Carr and Dritschel2010), is to ‘virtually’ extend the ambient density profile ${\it\rho}_{a}(z_{a})$ below $z_{a}=0$ by smoothly and rapidly matching the stratification at the boundary to a uniform-density region ${\it\rho}_{a}(-\infty )$ , only slightly greater than ${\it\rho}_{a}(0)$ , for $z-{\it\eta}<0$ . This results in a uniform-density core with approximately zero circulation. For example,

Incorporating this approximation into (3.2) allows smooth solutions for ${\it\eta}(z)$ throughout the domain. We employ this approximation with ${\it\delta}=0.01H$ only when $z-{\it\eta}<0$ somewhere in the domain. In practice, this occurs only where there is appreciable ambient stratification near $z=0$ (Lamb Reference Lamb2002).

3.6. Numerical solution method

The continuous model, (3.4) and (3.7), can be solved by a nested iterative approach as follows (from the inner- to the outermost levels). (i) Beginning with an initial value for ${\it\eta}_{z}(0)\equiv {\it\eta}_{o}^{\prime }$ , a proxy for the bore amplitude, find a value of $C_{b}$ that satisfies the eigenvalue problem (3.7) by a shooting method. That is, integrate the ordinary differential equation (ODE) with initial values ${\it\eta}(0)=0$ , ${\it\eta}_{z}(0)={\it\eta}_{o}^{\prime }$ , iterating on $C_{b}$ until the upper boundary condition, ${\it\eta}(H)=0$ , is satisfied. (ii) Iterate on ${\it\eta}_{o}^{\prime }$ until the bore amplitude, $h_{b}\equiv \int _{0}^{H}b(z)\,\text{d}z$ , is equal to the desired value. (iii) Iterate by varying the head loss constant, ${\it\Delta}_{o}$ , until the momentum conservation equation, (3.4), is satisfied.

4.1. Numerical method and set-up

To test the theory, numerical simulations were performed for internal bores released from a dam-break initial condition, with nominally two-layer stratification of varying interface thickness. We solve the two-dimensional Boussinesq Navier–Stokes equations,

where $p$ is the dimensionless pressure that remains after removing the hydrostatic part due to ${\it\rho}(H)$ , $b$ is the dimensionless buoyancy as defined in § 3.2, and $\hat{\boldsymbol{k}}$ is the unit vector in the positive $z$ direction. Velocities are scaled by $\sqrt{g^{\prime }H}$ , lengths by $H$ , time by $H/\sqrt{g^{\prime }H}$ , pressure by ${\it\rho}_{o}g^{\prime }H$ , and $\mathit{Re}$ is the Reynolds number. Model results presented below are scaled in the same manner.

The numerical model uses a finite-volume discretization with a second-order pressure projection method based on that of Bell & Marcus (Reference Bell and Marcus1992) and a Godunov-type advection scheme. This non-oscillatory finite-volume formulation is ideal for internal jumps with sharp gradients in density and velocity, and has been used to simulate gravity currents (White & Helfrich Reference White and Helfrich2008, Reference White and Helfrich2012) and nonlinear solitary waves (cf. Lamb Reference Lamb2002), correctly capturing phase speeds and exhibiting minimal energy loss over large distances.

Calculations were carried out in a rectangular domain spanning $-L/2\leqslant x\leqslant L/2$ and $0\leqslant z\leqslant H$ , where $L=32$ and $H=1$ , with a resolution ( $x\times z$ ) of $8192\times 256$ . Boundary conditions were free slip and a dam-break initial condition was used, with $\boldsymbol{u}=0$ and an initial buoyancy field

This results in a virtual dam of height $h_{d}$ , that transitions to an ambient (rescaled) density profile, ${\it\rho}_{a}(z)=(1/2)-\tanh [{\it\lambda}(z-h_{a})]/2$ with a nominal interface located at $h_{a}$ and a pycnocline thickness described by the parameter ${\it\lambda}$ (thickness ${\sim}1/{\it\lambda}$ ). Approximately 150 simulations were conducted spanning ${\it\lambda}=[4,8,12,24,64]$ , $r\equiv h_{a}/H=[0.1,0.2,0.3,0.4]$ and a range of dam heights between $h_{d}=[0,1]$ for each $({\it\lambda},r)$ combination. The Reynolds number for the simulations is $\mathit{Re}\equiv \sqrt{g^{\prime }H^{3}}/{\it\nu}=40\,000$ . While there is always a small amount of numerical dissipation, the resolved viscous dissipation, globally integrated, was typically 85 % of the total energy loss calculated from the residual of the kinetic energy budget.

4.2. Internal bore properties

In order to compare the theory with the numerical simulations, the bore speed $C_{b}$ and effective height $h_{b}$ were extracted from the results. The height is taken to be

which makes the hydrostatic pressure behind the bore equivalent to that of a sharp two-layer bore of the same height. This definition was also used by Borden et al. (Reference Borden, Meiburg and Constantinescu2012) for internal bores and by Marino, Thomas & Linden (Reference Marino, Thomas and Linden2005) for the equivalent thickness of a gravity current. The mean height, $\overline{h}_{b}$ , is the average value of $h_{b}(x)$ between $x=[0,x_{f}]$ , where $x_{f}$ is the front position, taken as the point at which the displacement is a quarter of the maximum, $h_{b}(x_{f})=[h_{b_{max}}-\tilde{h}_{a}]/4$ , where $\tilde{h}_{a}=h_{b}(L/2)$ . (Note that, because $\tilde{h}_{a}$ is the equivalent thickness, only in the sharp two-layer limit is it necessarily the case that $\tilde{h}_{a}=h_{a}$ , although in general they are very nearly the same.) The bore speed $C_{b}$ is calculated for each case by linear regression of $x_{f}$ against $t$ . From figure 2(b), it can be seen that the regression is highly linear; the uncertainty in $C_{b}$ is generally less than 0.1 % of its mean value.

Figure 2.

Internal bore properties from numerical simulations ( $h_{d}=0.4$ , $r=0.2$ , ${\it\lambda}=12$ shown as an example). (a) Bore height $h_{b}(x)$ and mean value $\overline{h}_{b}$ , taken between the initial dam location $x_{d}$ and the front position $x_{f}\ (t=20)$ . Bore height at the location of the first maximum is $\overline{h}_{b_{max}}$ . (b) Bore speed $C_{b}$ , obtained by linear regression of front position versus time.

Figure 3 shows the time evolution of an undular bore (a) and a turbulent bore (b) for $r=0.2$ and ${\it\lambda}=24$ . The undular bore, which evolves from an initial dam height $h_{d}=0.5$ , is transient and continually radiates solitary-like waves from the front (cf. Esler & Pearce Reference Esler and Pearce2011). The turbulent bore, with initial dam height $h_{d}=0.9$ , exhibits a smooth monotonic front with Kelvin–Helmholtz vortices and turbulent mixing behind. While there is a continuum between these states, depending on $h_{d}$ for a given $[{\it\lambda},r]$ , this dichotomy provides a reasonable means of classifying the simulations.

Figure 3.

Results from dam-break simulations shown at five different times (same for both panels): (a) undular bore ( $h_{d}=0.5$ , $r=0.2$ , ${\it\lambda}=24$ ); (b) turbulent bore ( $h_{d}=0.9$ , $r=0.2$ , ${\it\lambda}=24$ ). Greyscale (colour online) shows density.

Figures 4 and 5 illustrate the effect of interface thickness, i.e.  ${\it\lambda}$ , on undular and turbulent bores. For undular bores, $r=0.2$ , $h_{d}=0.5$ , a decrease in interface thickness affects the rate of internal wave radiation, producing an increasing number of crests with decreasing thickness (figure 4). For turbulent bores, perhaps better termed monotonic bores (figure 5), increasing interface thickness reduces the turbulent mixing, suppressing the Kelvin–Helmholtz instability entirely for the ${\it\lambda}=4$ case. The onset of a critical Richardson number, $\mathit{Ri}=0.25$ , may be a good measure to distinguish between the cases, as discussed in § 5.1.

Figure 4.

Undular bores with varying interface thickness, ${\it\lambda}$ , as shown: (a) density field; (b) velocity, density and Richardson number profiles from model, with ${\it\epsilon}\rightarrow -\infty$ (solid lines), and simulations (dashed lines).

Figure 5.

Turbulent (monotonic) bores with varying ${\it\lambda}$ , as shown. Details are as described in figure 4. Critical Richardson number ( $\mathit{Ri}=0.25$ ) shown by blue line.

5.1. Velocity and density profiles

The shear and stratification within the bore are important properties because together they influence the onset of shear instability and turbulent mixing. Profiles of mean velocity $U(z)$ , density $b(z)$ and gradient Richardson number $\mathit{Ri}=N^{2}/(\text{d}U/\text{d}z)^{2}$ are shown in figures 4 and 5 (right panels) along with model predictions, using ${\it\epsilon}\rightarrow -\infty$ . Mean profiles were calculated from the simulations by averaging over $x=[0,x_{f}]$ , and vertically shifting to align local profiles with the mean bore thickness. That is, we average $U(x,\tilde{z})$ and $b(x,\tilde{z})$ , where $\tilde{z}=z-(h_{b}(x)-\overline{h}_{b})$ . This shift means that undular motions do not contribute to pycnocline thickness, while still preserving local gradients, leaving shear, $N^{2}$ and $\mathit{Ri}$ unaltered. The instantaneous profiles were then averaged over time between $t=[10,T_{f}]$ , where $T_{f}$ is the final time for each simulation.

For the undular bores shown in figure 4, the model prediction based on ${\it\epsilon}\rightarrow -\infty$ and the simulations are almost indistinguishable. Here the Richardson number was always greater than unity, so it was not shown over the entire range. For turbulent (or monotonic) bores (figure 5), the agreement is still quite good. The minimum Richardson number is well predicted by the model, and shear instabilities are well correlated with minimum $\mathit{Ri}$ that approach $1/4$ . Note that the velocity and density profiles may be somewhat sensitive to the relative proportions of the turbulent to smooth frontal regions. These proportions may change to some degree as the bore propagates and both regions grow (see e.g. figure 3 b). However, the time averaging should to some degree account for these effects.

5.2. Bore speed

For comparison with the Navier–Stokes simulations, theoretical curves for $C_{b}(h_{b})$ were generated over a range of $({\it\lambda},r)$ values using the conjugate bore theory described in § 3. The effect of dissipation model was also studied by varying ${\it\epsilon}$ . Results are shown in figure 6 for $r=0.1$ (a–c), 0.2 (d–f), 0.3 (g–i) and 0.4 (j–l). The KRS two-layer model is shown for comparison. Each model curve is terminated at the energy-conserving conjugate state, ${\it\Delta}_{o}=0$ , which is independent of ${\it\epsilon}$ for a given $({\it\lambda},r)$ , since $D<0$ for larger amplitude.

Figure 6.

Bore speed versus amplitude. Comparison between model and theory for bores with varying ambient stratification. Symbols (with error bars) represent numerical results based on $C_{b}$ and $\overline{h}_{b}$ : ${\it\lambda}=4$ (○),  ${\it\lambda}=12$ (▫), ${\it\lambda}=24$ (▵). The energy-conserving conjugate state is denoted by $\ast$ . Solid symbols show the $h_{d}=0.5$ and $h_{d}=0.9$ cases for comparison with figures 4 and 5. The curves show the theoretical results for varying ${\it\epsilon}$ and ${\it\lambda}$ (see legends): KRS two-layer model; ${\it\epsilon}=1$ (continuous CBWS); ${\it\epsilon}=0$ (continuous VS); ${\it\epsilon}=-1$ (continuous KRS); ${\it\epsilon}\rightarrow -\infty$ (continuous maximal dissipation (MD) model).

In general, the agreement between the model and the simulations is excellent. The effect of the finite pycnocline is significant, in general decreasing $C_{b}$ for increasing pycnocline thickness. The model reproduces the correct small-amplitude limit, where the bore speed approaches the linear long-wave phase speed, $C_{b}\rightarrow c_{o}$ . For large-amplitude bores, all models (independent of ${\it\epsilon}$ ) tend to the conjugate-state limit, consistent with DJL nonlinear solitary-wave solutions, and demonstrating that the nonlinear ambient waveguide is structuring the bore speed. The two-layer theories, regardless of the model for dissipation or turbulent mixing, are incapable of reproducing the effect of ambient interface thickness on the waveguide.

For larger values of the ambient lower-layer thickness ( $r=0.3,\ 0.4$ ), the agreement is nearly exact up to the conjugate state, and, moreover, is nearly independent of ${\it\epsilon}$ . For $r=0.1,\ 0.2$ and intermediate values of the amplitude, the theory is more sensitive to the dissipation model, although the agreement with the simulations is still reasonable for most ${\it\epsilon}$ . As with two-layer bores, the continuous analogue of the KRS model ( ${\it\epsilon}=-1$ ) better captures the bore speed at large amplitude than the CBWS analogue ( ${\it\epsilon}=1$ ). However, the best agreement for all values of ( ${\it\lambda},r$ ) is found for the ${\it\epsilon}\rightarrow -\infty$ model.

The interface thickness has an even greater influence on the front speed than the dissipation model. With increasing thickness, the speed curves shift downwards, consistent with the ambient waveguide (e.g. decreasing long-wave speed $c_{o}$ and conjugate-state speed $C_{cs}$ ). Even a small but finite ambient pycnocline thickness has a significant effect on the bore speed relative to the two-layer limit. Note that Borden & Meiburg (Reference Borden and Meiburg2013b ) also found a decrease in bore speed for finite downstream pycnocline thickness, which their semi-empirical model captured relative to two-layer theories of KRS and CBWS and their own VS model. However, they do not address the role of finite pycnocline thickness on the ambient waveguide, which we show here to have a first-order effect on the bore speed.

5.3. Connection with the nonlinear internal waveguide and DJL theory

The transition from undular to monotonic bores with turbulent mixing can be viewed as a continuum of solitary waves of increasing amplitude, consistent with the DJL theory discussed in § 2.2 and also the dispersive two-layer theory of Esler & Pearce (Reference Esler and Pearce2011). There is a direct link between the bore speed, $C_{b}$ , and the intrinsic speed of nonlinear internal waves at the bore front. To illustrate this link, we have calculated the best-fitting solitary-wave solution to the DJL equation (2.1) and compared it to the leading wave at the bore front. DJL solutions can be characterized by their available potential energy, $\mathit{APE}=\iiint _{0}^{{\it\eta}(x,z)}b(z-{\it\eta}(x,z))-b(z-{\it\xi})\,\text{d}{\it\xi}\,\text{d}x\,\text{d}z$ (Stastna & Lamb Reference Stastna and Lamb2002). For comparison, we calculate from simulations the APE at the bore front, between the ambient and the point of the first maximum, $h_{b_{max}}$ . Because solutions to (2.1) are symmetric about the crest, we compare them with DJL solitary waves with an APE of twice this value. Results are shown in figure 7 for waves of increasing amplitude (increasing $h_{d}$ ) for $r=0.2,\ {\it\lambda}=8$ . Figure 8 shows the DJL predictions for solitary-wave speed as a function of amplitude compared against the simulation results using the bore speed and the maximum amplitude at the first crest. DJL solutions that develop trapped cores (see § 3.5) are extended in the core region using (3.14).

Figure 7.

Characteristics of the bore front compared with the DJL nonlinear solitary-wave theory. Greyscale plot (colour online) shows density field from simulations for bores of increasing amplitude. Solid line is the $b=0.5$ isopycnal from the DJL solution, obtained by matching APE for $x\geqslant h_{b_{max}}$ .

Figure 8.

Bore speed versus maximum amplitude, $h_{b_{max}}$ , compared with the DJL theory. DJL theory: solid lines. Numerical results: ${\it\lambda}=4$ (○), ${\it\lambda}=8$ (▫), ${\it\lambda}=24$ (▵). Solid squares: conjugate state. Solid circles: onset of trapped core solutions.

For both the front shape and the front speed, it is clear that the agreement between the simulations and DJL solutions is almost exact up to the conjugate state, the maximum amplitude for DJL solutions. This illustrates that the bore speed is intrinsically linked to the characteristics of the waves at the front, and the nonlinear DJL model is able to capture very accurately both the relationship $C_{b}(h_{b_{max}})$ and the shape of the front. Our hydraulic model for the bore speed can therefore be viewed as the function that connects $h_{b_{max}}$ to $\overline{h}_{b}$ . In this view, the upstream condition sets the bore speed, which then determines the characteristics of waves at the front, e.g. their amplitude and wavelength, and whether the bore is undular or monotonic/turbulent.

To further illustrate the connection between the bore model and the waves at the front, figure 9 shows the relationship between $h_{b_{max}}$ and $\overline{h}_{b}$ . These results can be compared against two limits. For a smooth monotonic bore, $h_{b_{max}}=\overline{h}_{b}$ , i.e. a 1:1 relationship, shown in figure 9. For undular bores, $h_{b_{max}}>\overline{h}_{b}$ , with a specific relationship that is determined by a complex balance between nonlinear and dispersive effects. However, in the weakly nonlinear Korteweg–de Vries (KdV) limit for a single-layer fluid, Whitham modulation theory predicts that the amplitude of the leading wave is twice the bore amplitude (Whitham Reference Whitham1974), and this 2:1 curve is also shown in figure 9. (Note that Esler & Pearce (Reference Esler and Pearce2011) extended the Whitham modulation theory to fully nonlinear weakly dispersive two-layer undular bores.) Most of the simulations for $r=0.1$ and $r=0.2$ result in $h_{b_{max}}>\overline{h}_{b}$ , i.e. undular bores. As the conjugate state is approached, bores become monotonic with $h_{b_{max}}\approx \overline{h}_{b}$ . The Whitham limit should apply for bores with a deep upper layer (small $r$ ), a thin interface (large ${\it\lambda}$ ) and small bore amplitude, so the most relevant comparison shown in figure 9 is for $r=0.1$ and ${\it\lambda}=24$ . Indeed, the small-amplitude solutions for this case initially follow the 2:1 line, but then diverge, eventually approaching the 1:1 line as the conjugate state is approached. The $r=0.3$ and $r=0.4$ cases each approximately follow the 1:1 line, suggesting nearly monotonic bores for all amplitudes.

Figure 9.

Maximum bore amplitude versus mean bore amplitude. Numerical results: ${\it\lambda}=4$ (○),  ${\it\lambda}=8$ (▫), ${\it\lambda}=24$ (▵). Solid symbols: conjugate state. Dashed line is $1:1$ relationship. Solid line represents the prediction of Whitham modulation theory in the single-layer KdV limit, $(h_{b_{max}}-h_{a})=2(h_{b}-h_{a})$ .

5.4. Energy dissipation

Results for the total energy dissipation, $D$ , from simulations and model predictions are compared in figure 10. The model results are calculated from (3.9) across a control volume applied between $x=[0,x_{f}]$ and subsequently averaged over time between $t=[10,T_{f}]$ . The model predictions are based on (3.10).

Figure 10.

Bore dissipation versus amplitude. Comparison between the model and theory for varying dissipation models (a,b) and interface thickness (c,d). Symbols represent numerical results based on measured $D$ and $\overline{h}_{b}$ : ${\it\lambda}=4$ (○), ${\it\lambda}=12$ (▫), ${\it\lambda}=24$ (▵). Solid symbols show the $h_{d}=0.5$ and $h_{d}=0.9$ cases for comparison with figures 4 and 5. In panels (a,b), curves show the theoretical results for varying ${\it\epsilon}$ , for (a)  $r=0.1$ and (b)  $r=0.2$ . In panels (c,d), curves show the theoretical results for varying ${\it\lambda}$ , for (c)  $r=0.1$ and (d)  $r=0.2$ .

Figure 10(a,b) compares the dissipation from simulations versus the model with varying ${\it\epsilon}$ . For all cases, the ${\it\epsilon}\rightarrow -\infty$ model produces the greatest net dissipation as compared with the CBWS and KRS analogues, despite the energy gain in the lower layer, as discussed in § 3.3. This explains the success of the ${\it\epsilon}\rightarrow -\infty$ model in predicting the bore speed, as discussed in § 5.2. The higher rate of dissipation reduces the bore speed (for solutions with ${\it\Delta}_{o}>0$ , up to the conjugate state), correcting the tendency of the other models to overpredict $C_{b}$ , which is also seen in the two-layer limit (see Borden et al. Reference Borden, Meiburg and Constantinescu2012).

For small-amplitude bores, the dissipation tracks the model solutions fairly well, falling either along the KRS or ${\it\epsilon}\rightarrow -\infty$ curves. However, for larger amplitudes, the dissipation is less than the model prediction, and becomes negative at a smaller amplitude than the model conjugate state (the $D=0$ solution). The negative dissipation signifies an energy gain over the region where the bore characteristics are measured. Note that this region encompasses only the rightward-propagating bore, and, as must be the case, when the leftward-travelling disturbance created by the initial dam break is also included in the calculation, there is a net energy loss over the entire volume.

Figure 10(c,d) shows the simulation results versus the ${\it\epsilon}\rightarrow -\infty$ model curves for a range of ${\it\lambda}$ . The effect of ${\it\lambda}$ is subtle. For $r=0.1$ the maximum dissipation occurs for the mid-range, ${\it\lambda}=12$ , in both the model and simulations for all amplitudes. For $r=0.2$ the model predicts that the maximum dissipation occurs at ${\it\lambda}=4$ , whereas the ${\it\lambda}=12$ simulation produces the maximum level of dissipation. Overall, the dissipation from the simulations is less than any of the model predictions for larger-amplitude bores. Nonetheless, the model captures the order of magnitude of the qualitative trends of dissipation with amplitude fairly well.

The net dissipation is due to a combination of turbulent mixing and wave drag, where the latter is due to correlations in pressure and velocity fluctuations that arise when averaging over the undular bore. Previous studies (e.g. Borden et al. Reference Borden, Meiburg and Constantinescu2012) have focused on turbulent mixing alone as the source of dissipation. However, it becomes clear from figure 10 that the large-amplitude turbulent bores, which are characterized by Kelvin–Helmholtz instabilities and mixing behind the front, have smaller (and in some cases negative) dissipation than the intermediate-amplitude bores, which often do not exhibit turbulent mixing. These results suggest that wave drag rather than turbulent mixing is the dominant mechanism of energy dissipation. This is consistent with the results of § 5.3, which showed that the internal wave characteristics at the bore front are of leading-order importance in determining the bore characteristics.

5.5. Velocity and head loss in isopycnal coordinates

The conjugate bore model, like the full DJL model, can be readily formulated in isopycnal coordinates. Moreover, the dissipation closure we have introduced casts the head loss as a function of buoyancy. It is therefore instructive to compare the model predictions and simulation results in isopycnal coordinates. Figure 11 shows profiles of velocity and Bernoulli head loss in this coordinate system.

Figure 11.

Profiles of velocity $U-C_{b}$ (first and third rows) and Bernoulli head loss (second and fourth rows) normalized as ${\it\Delta}/{\it\Delta}_{max}$ in isopycnal coordinates. Results correspond to the undular and turbulent bores from figures 4 and 5, with $r=0.2$ , $h_{d}=0.5$ and 0.9, and ${\it\lambda}=4$ , 12 and 24, as indicated. In all plots, solid symbols are the numerical results and lines are model predictions based on three different dissipation models for comparison: continuous VS model, ${\it\epsilon}=0$ (solid line); continuous KRS model, ${\it\epsilon}=-1$ (dashed line); and MD model, ${\it\epsilon}\rightarrow -\infty$ (dash-dotted line).

From the simulations, we calculate the local head loss in isopycnal coordinates from the Bernoulli function (3.1),

where $u(x,z)$ , $w(x,z)$ , $p(x,z)$ and $b(x,z)$ are the local fields, ${\it\eta}(x,z)=z-z_{a}(b)$ and $u$ is in the frame moving with $C_{b}$ . We calculate the average values of ${\it\Delta}$ between $[0,x_{f}]$ and between $t=[10,T_{f}]$ for each simulation. The velocity profiles are averaged from the simulations as described in § 5.1 to yield $U(b(\tilde{z}))$ .

The simulation results are compared with model predictions using three separate dissipation models: the continuous VS model, ${\it\epsilon}=0$ ; the continuous KRS model, ${\it\epsilon}=-1$ ; and the dissipation maximizing model, ${\it\epsilon}\rightarrow -\infty$ . For the velocity profiles, the model–simulation agreement is very good for all cases. The effect of the dissipation models is almost indistinguishable, despite the fact that the model choice has a demonstrable effect on $C_{b}$ (see figure 6).

The head loss, for both model and simulations, is normalized by the maximum value, ${\it\Delta}_{max}$ . This is instructive because, for each of the three models, ${\it\epsilon}=[0,-1,-\infty ]$ , the maximum loss occurs for $b=0$ , which corresponds to the upper boundary, $z=H$ . Based on this normalization, the differences between the VS model (uniform ${\it\Delta}$ ), the KRS model (energy conservation at the lower boundary, ${\it\Delta}(b=1)=0$ ) and the ${\it\epsilon}\rightarrow -\infty$ model (energy gain at the lower boundary, ${\it\Delta}(b=1)=-{\it\Delta}(b=0)$ ) are clear.

The undular bores from the simulations, shown in figure 11(b,d, f), show a consistent trend. For each, ${\it\Delta}$ is positive for $b=0\ (z=H)$ and decreases with increasing $b$ (decreasing $z$ ). For each, there is a region where ${\it\Delta}<0$ near $b=1$ , corresponding to an energy gain in the expanding region of the bore, near $z=0$ . There is also a region where ${\it\Delta}>0$ near $b=0$ , corresponding to an energy loss in the contracting region of the bore, near $z=1$ . The results suggest a clear energy transfer from the contracting to expanding regions for undular bores. None of the profiles is linear in $b$ , as the dissipation models postulate, although each is approximately linear for intermediate isopycnal values. Moreover, the ${\it\epsilon}\rightarrow -\infty$ model is arguably the closest match as compared with the KRS and VS models, with the best fit lying somewhere between ${\it\epsilon}=-1$ and ${\it\epsilon}\rightarrow -\infty$ .

The turbulent/monotonic bores from the simulations, shown in figure 11(h, j,l), show a less coherent trend. The ${\it\lambda}=4$ and 12 cases each have a distinct peak near an intermediate isopycnal value, more likely associated with mixing at the interface. The ${\it\lambda}=24$ case shows a structure that is more similar to the undular bores, despite the fact that it is a fully turbulent case, for which interfacial mixing is high. The profiles of ${\it\Delta}$ vary considerably with interface thickness and bore amplitude, and in general are nonlinear in $b$ . Nonetheless, the model for ${\it\epsilon}$ does not seem to have a significant influence on the velocity structure. Recall that turbulent bores have smaller net dissipation than undular bores, and that the dissipation decreases (see figure 10) and the bore speed is less sensitive to ${\it\epsilon}$ (figure 6) as the conjugate state is approached. It may be that the ${\it\epsilon}\rightarrow -\infty$ model gives a better prediction of $C_{b}$ because it better captures the energy transfer from the contracting to the expanding region for the higher-dissipation undular bores.