2.1 Computational approach, domain and set-up

The Navier–Stokes equations in an inertial frame, using the Boussinesq approximation for an inhomogeneous, incompressible fluid subject to gravity along the vertical direction $z$ are used to simulate the laboratory-scale system. The stability of the bolus structure in the transverse direction as it propagates up the slope (Venayagamoorthy & Fringer Reference Venayagamoorthy and Fringer2007) gives credence to using the two-dimensional model. However, we will verify that the coherent structure analysis of a three-dimensional simulation produces a similar bolus to the two-dimensional case in § 3.3. The Boussinesq approximation neglects effects of density variation except in the buoyancy term, and is appropriate for buoyancy-driven flows with weak relative density variations around a reference value $\unicode[STIX]{x1D70C}_{00}$. This approximation has been extensively used in ocean models and in two-layer density systems with small relative density change (Long Reference Long1965; Helfrich & Melville Reference Helfrich and Melville2006; Pedlosky Reference Pedlosky2013). The system of equations is

(2.1)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.2)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}u+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})u=-\frac{1}{\unicode[STIX]{x1D70C}_{00}}\unicode[STIX]{x2202}_{x}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D735}^{2}u, & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}w+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})w=-\frac{1}{\unicode[STIX]{x1D70C}_{00}}\unicode[STIX]{x2202}_{z}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D735}^{2}w-\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{00}}g, & \displaystyle\end{eqnarray}$$

(2.4)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}^{2}\unicode[STIX]{x1D70C}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}=(u,w)$ is the velocity field, $p$ is the pressure, $\unicode[STIX]{x1D70C}$ is the local density, $\unicode[STIX]{x1D70C}_{00}=1000~\text{kg}~\text{m}^{-3}$ is the reference density, $\unicode[STIX]{x1D708}=10^{-6}~\text{m}^{2}~\text{s}^{-1}$ is the kinematic viscosity, $\unicode[STIX]{x1D705}=1.4\times 10^{-7}~\text{m}^{2}~\text{s}^{-1}$ is the thermal diffusivity for sea water (Kunze Reference Kunze2003) and $g=9.8~\text{m}~\text{s}^{-2}$ is the gravitational acceleration. It is assumed in (2.4) that density diffusion in the system is driven by thermal diffusion.

The system of equations (2.1)–(2.4) is solved numerically for $u$, $w$, $p$ and $\unicode[STIX]{x1D70C}$ using CDP 2.4, an unstructured, finite-volume based, large eddy simulation code that implements a fractional-step time-marching scheme (Ham & Iaccarino Reference Ham and Iaccarino2004; Mahesh, Constantinescu & Moin Reference Mahesh, Constantinescu and Moin2004). All subgrid-scale modellings are turned off, and the code is modified to include the buoyancy term in (2.3) and solve (2.4) along with (2.1)–(2.3). This code has previously been used to simulate internal waves and has been validated with experiments (King, Zhang & Swinney Reference King, Zhang and Swinney2009; Dettner, Paoletti & Swinney Reference Dettner, Paoletti and Swinney2013; Lee et al. Reference Lee, Paoletti, Swinney and Morrison2014; Paoletti, Drake & Swinney Reference Paoletti, Drake and Swinney2014; Zhang & Swinney Reference Zhang and Swinney2014; Allshouse et al. Reference Allshouse, Lee, Morrison and Swinney2016; Lee et al. Reference Lee, Allshouse, Swinney and Morrison2018).

An illustration of the domain is presented in figure 2. The system dimensions correspond to an available system used for experimental studies (Allshouse et al. Reference Allshouse, Lee, Morrison and Swinney2016), with similar length scales as in previous experimental (Sutherland et al. Reference Sutherland, Barrett and Ivey2013; Moore et al. Reference Moore, Koseff and Hult2016) and numerical (Venayagamoorthy & Fringer Reference Venayagamoorthy and Fringer2007; Arthur & Fringer Reference Arthur and Fringer2014) bolus investigations. The primary reference frame $O_{xz}$ is positioned at the left-bottom corner of the domain. At the left end, a wavemaker generates an internal wave that propagates towards the slope. A constant-slope topography is positioned at the right end of the system. The topography, with constant slope $s$, is positioned such that the mid-depth line ($z=H/2$) intersects the slope at a distance $L=2.87~\text{m}$ from the left boundary. This point on the slope is used as the origin for the rotated reference frame $\widehat{O}_{\widehat{x}\widehat{z}}$, where $\widehat{x}$ and $\widehat{z}$ represent the coordinates along and normal to the slope, respectively. Because the location of $\widehat{O}$ is held constant, there is a minimum slope of $s=0.070$ below which the topography would reach the inlet. For there to be a constant depth development region, the minimum slope studied here is 0.105.

Figure 2. Schematic diagram of the computational domain, dimensions, mesh resolution, density profile (in red) and fundamental mode forcing velocity profile (in blue). The distance from the inlet boundary on the left to the midpoint of the slope $L=2.87~\text{m}$ and the domain height $H=0.4~\text{m}$ are constant in all simulations, while $\unicode[STIX]{x1D6FF}$ and $s$ are varied in the parametric study. A sample background density profile $\unicode[STIX]{x1D70C}_{0}(z)$ (red) and inlet velocity profile $U(z)$ (blue) are presented corresponding to the sample simulation ($\unicode[STIX]{x1D6FF}=0.2~\text{m}$, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=20~\text{kg}~\text{m}^{-3}$, $s=0.176$, presented in § 2.2). The mesh resolution of the sample simulation is presented in a grey logarithmic scale with the highest resolution $\unicode[STIX]{x1D6E5}_{x}=2\times 10^{-4}~\text{m}$ in the breaking region (black) and the lowest resolution $3\times 10^{-3}~\text{m}$ in the evanescent region (light grey).

In order to efficiently obtain the necessary accuracy for a direct simulation, unstructured triangular multi-block meshes with different resolution zones were adopted. The variable mesh resolution is illustrated in figure 2. While the mesh resolution is lower in weakly stratified regions, higher resolutions are used in regions of internal wave propagation ($z\in [0.1,0.3]$  m), and the highest resolution is used in the overturning and breaking region on the slope. For the representative mesh in figure 2, the mesh size, defined as the local averaged distance between cell centres, varies from $\unicode[STIX]{x1D6E5}_{x}=3\times 10^{-3}~\text{m}$ in the sparse regions, to $10^{-3}~\text{m}$ at the internal wave propagation zone, to a maximum resolution of $2\times 10^{-4}~\text{m}$. The meshes used for the simulations presented in this work contain between 3.8 and 5.5 million cells. Approximately 70 % of these cells are in the breaking zone (corresponding to the darkest zone in figure 2), where the complex velocity field needs to be accurately resolved. Simulations with the narrowest pycnocline thickness required a higher resolution mesh to provide converged results, with a total of 12.6 million cells and the smallest cell size in the breaking zone reduced to $5\times 10^{-5}~\text{m}$.

Convergence studies to define the required spatial and temporal resolutions have been conducted varying both the mesh resolution in the breaking zone and the time step used in the simulations. This study resulted in the use of a constant time step of $5\times 10^{-4}~\text{s}$ for all simulations. With velocity magnitudes below $0.05~\text{m}~\text{s}^{-1}$ and the smallest mesh size of $5\times 10^{-5}~\text{m}$, the maximum Courant number is 0.5. The sea water values for the kinematic viscosity and thermal diffusivity (Kunze Reference Kunze2003) correspond to a Prandtl number of $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}=7.14$. The Reynolds number in the bolus forming region can be estimated using a characteristic length $L_{c}=0.02~\text{m}$ and a characteristic speed $U_{c}=0.03~\text{m}~\text{s}^{-1}$ based on the bolus average size and speed, respectively, which corresponds to $Re=U_{c}L_{c}/\unicode[STIX]{x1D708}=600$. Therefore, the Kolmogorov length scale in the breaking region is $\unicode[STIX]{x1D702}=L_{c}Re^{-3/4}\approx 1.65\times 10^{-4}~\text{m}$. With a default mesh resolution of $\unicode[STIX]{x1D6E5}_{x}=2\times 10^{-4}~\text{m}$ in the breaking zone, of the same order as the microscale $\unicode[STIX]{x1D702}_{k}$, and $5\times 10^{-5}~\text{m}$ for the sharpest pycnocline thickness case, the breaking region is sufficiently resolved to capture the bolus dynamics.

The background density stratification is a continuous, decreasing function of $z$. For all simulations, a hyperbolic tangent profile of variable pycnocline thickness is used to provide a smooth transition between the densities at the top and bottom, representative of what is observed in the oceans (figure 1). The stratification is modelled as

(2.5)$$\begin{eqnarray}\unicode[STIX]{x1D70C}_{0}(z)=\unicode[STIX]{x1D70C}_{H/2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\tanh \left[\frac{2(z-H/2)}{\unicode[STIX]{x1D6FF}}\tanh ^{-1}(0.95)\right],\end{eqnarray}$$

for $z\in [0,H],$ where $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\approx \unicode[STIX]{x1D70C}_{0}(0)-\unicode[STIX]{x1D70C}_{0}(H)$ is the density change and $\unicode[STIX]{x1D70C}_{H/2}=\unicode[STIX]{x1D70C}_{0}(H/2)$ is the mid-depth unperturbed density value. The parameter $\unicode[STIX]{x1D6FF}$ is the pycnocline thickness for the hyperbolic tangent profile and corresponds to the transition height in which the density varies by 95 % of $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$, as illustrated in figure 3(a). The hyperbolic tangent profile shape has been previously used to model sharp density stratifications in the limiting case of an almost two-layer density fluid (corresponding to small $\unicode[STIX]{x1D6FF}$) because of its stability properties compared to a discontinuous profile (Thorpe Reference Thorpe1971; Fringer & Street Reference Fringer and Street2003; Troy & Koseff Reference Troy and Koseff2005; Arthur et al. Reference Arthur, Koseff and Fringer2017). Because the pycnocline is located at the centre of the domain, two-layer theory does not predict mode-1 internal solitary wave propagation (Long Reference Long1956). Unlike Arthur & Fringer (Reference Arthur and Fringer2014), we will not initialize the wave with a form similar to the solitary wave solution, and instead force the system with a vertical mode-1 profile.

Figure 3. (a) Background density stratification for a density change $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=20~\text{kg}~\text{m}^{-3}$ and pycnocline thicknesses $\unicode[STIX]{x1D6FF}=0.025$ (blue), 0.2 (red) and 0.4 m (green), as well as $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=10~\text{kg}~\text{m}^{-3}$ and $\unicode[STIX]{x1D6FF}=0.2~\text{m}$ (yellow). (b) The normalized mode-1 velocity profile amplitude $U(z)$ for each corresponding stratification, with frequency $\unicode[STIX]{x1D714}=0.628~\text{rad}~\text{s}^{-1}$. The red horizontal dashed lines indicate the pycnocline thickness for $\unicode[STIX]{x1D6FF}=0.2~\text{m}$, and the yellow vertical dashed lines the density change for $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=10~\text{kg}~\text{m}^{-3}$.

At $t=0$, the unperturbed system is at rest ($u=w=0$), the density field is $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}(z)$ and the hydrostatic pressure field is $p=p_{0}(z)$. As illustrated on the left boundary of the domain in figure 2, the system is perturbed from its quiescent state by forcing the horizontal velocity $u$ of the vertical mode-1 plane wave of frequency $\unicode[STIX]{x1D714}=0.628~\text{rad}~\text{s}^{-1}$ (10 s-period). The forcing mechanism reproduces numerically what would be obtained experimentally by an oscillating plate wavemaker (Mercier et al. Reference Mercier, Martinand, Mathur, Gostiaux, Peacock and Dauxois2010). The mode-1 wave is particularly interesting as it is known that internal wave generation mechanisms create waves with the majority of the energy in the lowest modes, with the first mode being the one with the largest wavelength and lowest shearing stress, and therefore the least affected by viscous dissipation (Gerkema & Zimmerman Reference Gerkema and Zimmerman2008). The finite-difference approach for determining the vertical mode-1 profile of frequency $\unicode[STIX]{x1D714}$ (and subsequent modes) for an arbitrary stratification is presented in the supplementary material available at https://doi.org/10.1017/jfm.2019.993.

For $\unicode[STIX]{x1D70C}_{0}(z)$ in the form (2.5), the vertical mode $U(z)$ is a function of $\unicode[STIX]{x1D714}$, $H$, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D6FF}$. Figure 3(b) presents how the fundamental mode profiles $U(z)$ differ for density profiles $\unicode[STIX]{x1D70C}_{0}(z)$ varying the value of $\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$. Note that smoother density profiles correspond to smoother velocity profiles. Thinner pycnoclines are associated with a stronger velocity shear, $\text{d}U/\text{d}z$, around $z=H/2$. The left boundary condition for $u$ is given by

(2.6)$$\begin{eqnarray}u(x=0,z,t)=\left\{\begin{array}{@{}ll@{}}a_{f}\,U(z)\sin (\unicode[STIX]{x1D714}t)\quad & \text{for }t\in [0,2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}],\\ 0\quad & \text{for }t>2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714},\end{array}\right.\end{eqnarray}$$

where $U(z)$ is normalized to unit maximum magnitude and $a_{f}$ prescribes the amplitude of the forcing. The system is forced for a single period, because we are only interested in the shoaling of the first, leading wave. After one period, the velocity at the left boundary is set to zero. Only the horizontal velocity is modified from the background value for the forcing. Because the system is forced with a mode-1 wave, we are not able to independently modify properties of the wave such as wave speed, amplitude or wavelength. Additionally, $a_{f}$ is large enough to cause the propagating wave to be nonlinear, further complicating the relationship between these properties and the wave forcing parameters.

Changing the stratification profile modifies the shape of the velocity profile $U(z)$ forced on the left boundary (as illustrated in figure 3), but it does not prescribe $a_{f}$. To compare similar waves while varying parameters, we want the kinetic energy of the breaking waves to be held constant, and the forcing amplitude $a_{f}$ is thus set by the amount of energy present in the resulting breaking wave. The wave front kinetic energy through the vertical transect at the horizontal position $x$ and time $t$ is quantified as

(2.7)$$\begin{eqnarray}E_{k}(x,t)=\int _{0}^{H}\frac{1}{2}\unicode[STIX]{x1D70C}(x,z,t)[u^{2}(x,z,t)+w^{2}(x,z,t)]\,\text{d}z,\end{eqnarray}$$

which has units of energy per unit area. As the wave travels through the domain, it experiences dissipation and dispersion, changing shape, amplitude and speed even before it arrives at the start of the slope. Such changes depend on the stratification profile and on the amplitude of the forcing, in such a way that injecting waves at the inlet with the same maximum front kinetic energy $E_{k}(0,\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714})$ was found not to be equivalent to obtaining waves of equal front kinetic energy at the breaking location. To get the same kinetic energy at the breaking location, preliminary simulations are run in a constant depth channel to determine the forcing velocity amplitude $a_{f}$. This process guarantees that the generated wave fronts have the same instantaneous kinetic energy at the breaking point $x=L$, at the time $t_{c}(L)$ when the leading wave crest is at $x=L$. The value used for the kinetic energy as defined in (2.7) is $E_{k}(L,t_{c}(L))=E_{k,0}=6.734\times 10^{-3}~\text{J}~\text{m}^{-2}$ for all cases.

The boundary conditions for the top, bottom and sloping boundaries are no slip for velocity, $u=w=0$, and no flux for density. Imposing a no-slip condition at the top boundary, which corresponds to the water surface, does not significantly impact the dynamics because internal wave perturbations are strongest around mid-depth and decrease exponentially towards top and bottom, so that induced vertical velocities rapidly decay. To ensure numerical convergence and complement the inlet boundary condition on the left, a small ($5\times 10^{-3}~\text{m}$-high) outlet is added to the top-right of the domain, guaranteeing zero net flux within the domain, as expected from (2.1). Studies were performed to ensure that the presence, size and position of the outlet did not influence the dynamics in the breaking region.

2.2 Density perturbation field and bolus formation on slope

Before analysing the simulation results from a Lagrangian perspective, we present a typical bolus simulation to establish a basic intuition about the life cycle of the bolus. The density perturbation field, $\unicode[STIX]{x1D70C}^{\prime }(x,z,t)=\unicode[STIX]{x1D70C}(x,z,t)-\unicode[STIX]{x1D70C}_{0}(z)$, provides an Eulerian description of the wave propagation, breaking, bolus formation and propagation on slope. This field primarily corresponds to density fluctuations resulting from internal wave propagation, and because the bolus corresponds to higher density fluid moving up the slope, it will be highlighted by this field. The sample simulation parameters are $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=20~\text{kg}~\text{m}^{-3}$, $\unicode[STIX]{x1D6FF}=0.2~\text{m}$ and $s=0.176$. The mode-1 forcing profile for this case has a velocity amplitude of $a_{f}=0.0126~\text{m}~\text{s}^{-1}$, corresponding to oscillating plate displacements up to $0.02~\text{m}$ at frequency $\unicode[STIX]{x1D714}=0.628~\text{rad}~\text{s}^{-1}$. This forcing amplitude produces a front kinetic energy at the breaking site of $E_{k}(L,t_{c}(L))=E_{k,0}$.

Figure 4. (a) Time evolution of the mid-depth isopycnal, $\unicode[STIX]{x1D70C}_{H/2}=1010~\text{kg}~\text{m}^{-3}$, with the wave amplitude scale bar indicated in red. The dashed lines represent the unperturbed isopycnal height. (b) Instantaneous density perturbation field $\unicode[STIX]{x1D70C}^{\prime }$ of the sample simulation ($\unicode[STIX]{x1D6FF}=0.2~\text{m}$) at $t=55~\text{s}$ for the full domain and (c) breaking region. Positive values (red) represent fluid displaced upward and negative values (blue) represent fluid displaced downward. The solid grey lines in (b) represent isopycnals. The isopycnal $\unicode[STIX]{x1D70C}=1005~\text{kg}~\text{m}^{-3}$ is drawn in black and used to highlight the bolus front boundary in (c,d). The breaking region presented in (c) corresponds to the region surrounded by the thick black box in (b). (d) Velocity field corresponding to the region surrounded by the thick black box in (c), which contains the bolus front. (Video of the evolution of (b,c) is available in the supplementary material.)

The wave propagation and the density perturbation field are presented in figure 4. Figure 4(a) presents the time evolution of the mid-depth isopycnal $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{H/2}$, from $t=10~\text{s}$, after a full period of forcing, to $t=38~\text{s}$, when the isopycnal impinges on the topographic slope at $x=L$. The resulting wave is asymmetric, and its shape slowly changes as it propagates toward the slope. Note that despite exciting the system for a single period, several higher mode waves of decreasing amplitude are produced as a result of the initial forcing and persist even after the left boundary conditions are set to zero perturbation (for $t>2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}=10~\text{s}$). However, these waves move slower than the mode-1 wave and are rapidly left behind and do not affect the bolus dynamics.

The density perturbation field for the sample simulation at time $t=55~\text{s}$ is presented in figure 4(b,c). In this figure, white indicates regions of negligible density perturbation, red indicates positive density perturbation (i.e. denser fluid elements perturbed upwards from their equilibrium position) and blue indicates negative density perturbation (i.e. less dense fluid elements perturbed downwards). Figure 4(b) presents the full fluid domain: the generated waves (with alternating crests and troughs indicated in red and blue, respectively) have propagated to the right from the inlet and reached the slope. Small asymmetries in the generated waves (the maximum perturbations are not perfectly centred at $z=H/2$) demonstrate that the velocity amplitude is large enough to produce weakly nonlinear waves, even before they reach the sloping region. The solid grey lines in figure 4 represent isopycnals, and the isopycnal $\unicode[STIX]{x1D70C}=1005~\text{kg}~\text{m}^{-3}$ (black solid line) highlights the bolus front. The negligible density perturbation above and below the pycnocline demonstrates how the density perturbation amplitudes decrease exponentially in the vertical direction within the weakly stratified, evanescent regions.

The bolus resulting from the internal wave breaking consists of dense fluid moving upslope and is highlighted in figure 4(c) as a positive (red) density perturbation. The results are presented in the rotated frame $\widehat{O}_{\widehat{x}\widehat{z}}$, and the domain in figure 4(c) corresponds to the bold, black rectangle (0.9 m-long by 0.02 m-tall) in figure 4(b). The Eulerian results for the bolus propagation are qualitatively similar to those numerically produced by Venayagamoorthy & Fringer (Reference Venayagamoorthy and Fringer2007). The internal wave breaking forms the bolus, which corresponds to a propagating counter-clockwise vortex, as illustrated in figure 4(d) by the velocity field around the bolus front. The bolus dynamics as it enters the zone of less dense fluid is well described by the theory of gravity currents propagating in a stratified fluid (Benjamin Reference Benjamin1968; Maxworthy et al. Reference Maxworthy, Leilich, Simpson and Meiburg2002; White & Helfrich Reference White and Helfrich2008). The propagation of the dense front entering a zone of lower density fluid is subject to the induced shear and large difference in density between the bolus and the surrounding fluid, which results in vortex shedding and mixing that eventually stops the bolus from moving upslope. Towards the end of the bolus propagation, the bolus dynamics is impacted by other boluses produced by the following waves reaching the slope. These secondary boluses will not be addressed in this manuscript. We focus on the dynamics of the first bolus entering a quiescent region with no previous flow that will impact the bolus propagation.

While the density perturbation field and the velocity field provide good intuition for the bolus dynamics, a non-zero value of $\unicode[STIX]{x1D70C}^{\prime }$ at a given point only means that there is no local perturbation at that instant. This approach does not track fluid elements moving with the bolus, nor does this approach provide accurate information about how fluid is transported as the breaking happens. This is most acutely highlighted by the fact that the initial internal wave does not transport fluid from the inlet to the sloping boundary. At some point during the breaking process, the internal wave does begin to transport fluid (see the supplementary material for a demonstration video), and identifying when, how, and what fluid is transported is the purpose of the Lagrangian analysis.

3.1 Lagrangian characterization of the bolus

Lagrangian coherent structures provide a robust means for identifying the key underlying transport features of a given flow (Haller Reference Haller2002). One type of structure that can be detected is materially coherent vortices (Abernathey & Haller Reference Abernathey and Haller2018). These features are unique in that as they move through the domain, the fluid inside the region deforms but does not significantly mix with the fluid outside the perimeter. Identifying this type of feature in the breaking region will determine the fluid that is being advected by the bolus. Materially coherent features can be identified using Cauchy–Green-based metrics (Haller & Beron-Vera Reference Haller and Beron-Vera2013), transport operator methods (Froyland, Santitissadeekorn & Monahan Reference Froyland, Santitissadeekorn and Monahan2010), braid-based methods (Allshouse & Thiffeault Reference Allshouse and Thiffeault2012) and graph Laplacian methods (Froyland & Junge Reference Froyland and Junge2018). We use the clustering-based approach based on Hadjighasem et al. (Reference Hadjighasem, Karrasch, Teramoto and Haller2016) with minor modifications.

Central to coherent structure detection is the analysis of trajectories, which are computed from the velocity field $\boldsymbol{u}=(u,w)$ obtained by numerically solving (2.1)–(2.4). The velocity $\widehat{\boldsymbol{u}}$ in the rotated frame is used to advect massless fluid elements referred to as passive tracers, which have no inertia and move according to

(3.1)$$\begin{eqnarray}\displaystyle \frac{\text{d}\widehat{\boldsymbol{x}}_{i}}{\text{d}t}=\widehat{\boldsymbol{u}}(\widehat{\boldsymbol{x}}_{i}(t),t), & & \displaystyle\end{eqnarray}$$

where $\widehat{\boldsymbol{x}}_{i}(t)$ is the trajectory of tracer $i$ with initial position $\widehat{\boldsymbol{x}}_{i}(t_{0})$. Numerical integration of (3.1) is performed using a fourth-order Runge–Kutta method with a time step of $5\times 10^{-3}~\text{s}$. The breaking region is initially covered by a rectangular grid of passive tracers with along-slope dimensions that vary based on the slope and injected energy (the grids used for each simulation are in the supplementary material). The tracer spacing is $\unicode[STIX]{x0394}\widehat{x}_{p}=5\times 10^{-3}~\text{m}$ and $\unicode[STIX]{x0394}\widehat{z}_{p}=2\times 10^{-4}~\text{m}$, resulting in $n=18\,600$ to 40 000 tracers. The advection time window starts at time $t_{0}$, which is just before the shoaling wave arrives at the breaking region, and ends at time $t_{f}$, which corresponds to when the initial Eulerian density perturbation has stopped propagating up the slope, an upper bound for when the bolus stopped moving upslope.

Having computed the $n$ trajectories for the interval $[t_{0},t_{f}]$, we perform the coherent structure analysis based on the spectral clustering method of Hadjighasem et al. (Reference Hadjighasem, Karrasch, Teramoto and Haller2016). The clustering problem attempts to partition the domain into clusters such that trajectories within the same cluster are similar and trajectories from different clusters are dissimilar. The spectral analysis relies on the construction of a similarity graph that quantifies the pairwise similarity of trajectories. The similarity metric $w_{ij}$ is the inverse of the time-averaged distance between these trajectories. This calculation is performed for all pairs of trajectories, but only values of $w_{ij}$ greater than a user defined value are retained in order to sparsify the similarity matrix, $\unicode[STIX]{x1D652}$. Here, only values of $w_{ij}>1/r^{\ast }$, with $r^{\ast }=0.07~\text{m}$, are retained, resulting in an approximately 90 % sparse matrix for the sample simulation discussed in § 2.2. The clustering results were verified to be consistent for different choices of $r^{\ast }$ in a neighbourhood of this prescribed value. Based on the similarity matrix $\unicode[STIX]{x1D652}$, the diagonal degree matrix $\unicode[STIX]{x1D63F}$ is produced, where each diagonal element is equal to the sum of the elements in the corresponding row of $\unicode[STIX]{x1D652}$. From the sparse similarity matrix and degree matrix, the unnormalized graph Laplacian $\unicode[STIX]{x1D647}=\unicode[STIX]{x1D63F}-\unicode[STIX]{x1D652}$ is computed.

To create the clustering partition, we must identify characteristics of the set of trajectories. To characterize the trajectories, the next step of the algorithm is to compute the first generalized eigenvectors of the generalized eigenproblem

(3.2)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D647}\boldsymbol{q}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D63F}\boldsymbol{q}. & & \displaystyle\end{eqnarray}$$

It has been shown that the generalized eigenvalues of (3.2) satisfy $0=\unicode[STIX]{x1D706}_{1}\leqslant \cdots \leqslant \unicode[STIX]{x1D706}_{n}$, and the normalized dominant eigenvectors $\boldsymbol{q}_{1},\boldsymbol{q}_{2},\ldots ,\boldsymbol{q}_{n}$ differentiate properties in the graph and facilitate the clustering process (von Luxburg Reference von Luxburg2007). The dominant eigenvectors corresponding to the smallest eigenvalues reveal the most important characteristics of the flow. Each trajectory is characterized by a value within each eigenvector, and this characterization is used to group similar trajectories together. While all of the eigenvectors provide information, the dominant ones highlight the most important patterns.

The final step of the spectral clustering algorithm is to construct the eigenvector matrix $\boldsymbol{Q}\in \mathbb{R}^{n\times k}$ containing the $k$ dominant eigenvectors $\boldsymbol{q}_{1},\ldots ,\boldsymbol{q}_{k}$ as columns. While Hadjighasem et al. (Reference Hadjighasem, Karrasch, Teramoto and Haller2016) use the eigengap heuristic to determine $k$, this heuristic is not ideal for our application given the small amount of mixing outside of the bolus, so we adopted here a different heuristic based on the form of the dominant eigenvectors to select $k$. $\boldsymbol{q}_{k}$ is the first eigenvector highlighting a small fraction of tracers that entrain in a vortex and move upslope. Let $\boldsymbol{y}_{i}\in \mathbb{R}^{k}$ be the characterization vector corresponding to the $i$th row of $\boldsymbol{Q}$, which contains condensed differentiating information for trajectory $i$ (note that $k\ll n$). The characterization vectors $(\boldsymbol{y}_{i})_{1\leqslant i\leqslant n}$ are clustered with a $K$-means algorithm, assigning each vector $\boldsymbol{y}_{i}$ and the corresponding trajectory to a cluster. We use $k+1$ clusters to partition the domain, where an extra cluster is added to account for the incoherent cluster, as suggested by Hadjighasem et al. (Reference Hadjighasem, Karrasch, Teramoto and Haller2016).

Figure 5. (a) Sample simulation clustering result for elements in the initial uniform grid. Seven clusters have been identified, with the Lagrangian bolus cluster represented in dark blue. (b) The time evolution for the bolus cluster from $t_{0}=19.25~\text{s}$ to $t_{f}=65~\text{s}$, which is the entire bolus lifespan. (Video of the evolution of the figure is available in the supplementary material.)

Figure 5(a) presents the initial position of the partitions assigned to each one of the seven clusters identified for the sample simulation discussed in § 2.2. The dark blue cluster is identified by the method as the objective bolus: the one cluster composed of tracers that eventually entrain in a propagating vortex and move upslope. The bolus cluster also corresponds to the cluster of maximum displacement of the centre of volume in the $\widehat{x}$ direction. Figure 5(b) presents how the bolus cluster propagates in time, and we note that tracers maintain their cluster membership throughout the duration of the time interval. The bolus consists of tracers that are initially spread horizontally along the slope. As the wave arrives, those tracers get lifted, are trapped in the compact vortex, move along with the vortex for approximately 0.4 m up the slope, the front tracers stagnate and the trailing tracers start to recede down slope. All the elements identified as part of the bolus end up trapped inside of the vortex in intermediate times. While not the case for this example, it is possible that no bolus is detected by the algorithm, indicating that the shoaling internal wave does not result in effective transport.

3.2 Bolus transport properties

The objective quantification of the bolus makes it possible to measure bolus transport properties such as size, shape, position of centre of volume and velocity. Important bolus properties are here defined, applied to the sample simulation and will be used in the parametric studies in § 4. Position and velocity of the Lagrangian bolus propagating up the slope are also compared to the unbroken portion of the shoaling internal wave.

Let $I_{b}$ be the set of tracers identified as part of the bolus, so that $|I_{b}|=n_{b}$ is the number of passive tracers identified as the bolus. The tracer positions $\{\boldsymbol{x}_{i}\}$ are known in the uniformly sampled time interval $[t_{0},t_{f}]$. For each time instance $t$, averaging the position of the tracers gives the position of the bolus centre of volume

(3.3)$$\begin{eqnarray}\boldsymbol{x}_{CV}(t)=\frac{1}{n_{b}}\mathop{\sum }_{i\in I_{b}}\boldsymbol{x}_{i}(t),\end{eqnarray}$$

which we will use to quantify the bolus trajectory for $t\in [t_{0},t_{f}]$. The positions of the horizontally foremost and trailing tracers inside the bolus, $\boldsymbol{x}_{+}(t)$ and $\boldsymbol{x}_{-}(t)$ respectively, are also tracked. These positions will correspond to different tracers as time evolves and overturning inside the bolus takes place. These properties allow us to track the bolus centre of volume and its horizontal extent from a Lagrangian perspective.

To understand the relationship between the bolus and the shoaling wave that produces it, the wave trough and crest are also tracked. Because the unbroken components of the shoaling do not transport material, the displacement of isopycnals is tracked to locate the wave crest and trough. While the propagating wave amplitude is maximized at the centre of the pycnocline, other isopycnals throughout the water column are also deformed accordingly with these deformations vertically aligned. To track the unbroken shoaling wave, we track an isopycnal that is well above the breaking region, which makes it possible to track the shoaling wave propagation smoothly even after the bolus is generated. The first propagating minimum and maximum of the isopycnal correspond to the leading wave trough and crest horizontal position. Tracking the position of the isopycnal minimum and maximum provide a full description of the horizontal position of trough and crest as a function of time.

Figure 6. Kinematic comparison between the bolus and a co-propagating non-breaking isopycnal above the breaking region for $t\in [t_{0},t_{f}]$. (a) Bolus and amplified isopycnal for three time instances. The wave vertical positions and amplitudes are modified for illustration purposes. (b) The horizontal position as a function of time and (c) the horizontal speed as a function of horizontal position for the leading wave trough (dashed blue), the leading wave crest (dashed red), the bolus centre of volume (bold black), bolus trailing point (thin black) and foremost point (bold dark red). The isopycnal used to track the crest and trough in this plot was $\unicode[STIX]{x1D70C}=1000.05~\text{kg}~\text{m}^{-3}$, which has a neutrally buoyant height at $z=0.3635~\text{m}$. The tick marks in (b) correspond to the times for which the bolus and isopycnals are plotted in (a).

By tracking the extent of the bolus and the wave crest and trough, we can unveil how the location of the bolus relative to the wave evolves with time. Bolus and wave kinematics are compared in figure 6. A magnified version of the tracked isopycnal is presented in figure 6(a) at three different time instances, together with the corresponding boluses, demonstrating that the front of the bolus is always between the wave crest and trough. Figure 6(b) presents the $(x,t)$ curves for the horizontal extent of the bolus and the shoaling wave. For the bolus, the trajectories presented are those of the centre of volume, foremost and trailing point, with the horizontal extent of the bolus in grey. For the leading wave, the trough and crest trajectories are presented. As the wave trough approaches the slope, the trailing point of the bolus initially recedes while the front remains unaffected. The wave crest then starts to approach, pushing the rear of the bolus forward and ultimately creates a vortex, which entrains all of the bolus tracers. Finally, as depicted in figure 6(a) for $t=52~\text{s}$, the bolus extent begins to increase, as some of the tracers are shed from the vortex, while the foremost tracers continue moving upslope with the decelerating wave crest. While the leading edge tracers remain ahead of the wave crest, the ejection of tracers from the vortex and their recession downslope causes the centre of volume to be left behind. Eventually the trailing points begin to move back up slope, at $t\approx 62~\text{s}$, as the next shoaling wave arrives.

The horizontal speeds of the bolus centre of volume, foremost point and the leading crest are presented in figure 6(c). While the centre of volume moves slower than the crest, the correlated motion of the wave crest and the leading points is evident in the similar speed profiles with both undergoing the same deceleration as they move upslope where the water column narrows. Such a strong correlation between the Lagrangian bolus and the wave crest is not necessarily a result one would expect: the crest position is directly extracted from the density field, well above the breaking region, and the bolus is quantified by the movement of passive tracers identified as a coherent structure with the spectral clustering method. Nevertheless, the kinematics of both structures seem to consistently match, not only for the sample case presented, but for all parameter sweeps, which is an indicator that shoaling internal wave signatures close to the surface may possibly be used to estimate bolus location and speed. Given the consistency of the correlation between the crest and the bolus leading edge throughout the different cases analysed, we will not demonstrate it for each parameter sweep.

Instead, for a direct comparison between simulations with different parameters, we quantify how far upslope material is transported by using the trajectory of the bolus centre of volume in the rotated reference frame, $\widehat{x}_{CV}(t)$. The maximum displacement upslope relative to the initial position, $D_{b}$, is

(3.4)$$\begin{eqnarray}D_{b}=\max _{t\in [t_{0},t_{f}]}\{\widehat{x}_{CV}(t)-\widehat{x}_{CV}(t_{0})\}.\end{eqnarray}$$

For the sample case in figure 5(b), $D_{b}$ equals 0.402 m and is achieved at $t=56~\text{s}$. One could, alternatively, define this distance based on the position of the bolus foremost point, and that definition would provide a similar displacement upslope (0.363 m for the sample case). We consider the centre of volume to provide a better description of the Lagrangian bolus dynamics as it accounts for the complete coherent structure describing material trapped in the bolus.

Another quantity of interest is the amount of material being transported, which we quantify by the approximate area of the bolus. The bolus size, $S_{b}$, is estimated using the number $n_{b}$ of tracers in the bolus cluster,

(3.5)$$\begin{eqnarray}S_{b}=n_{b}\,\unicode[STIX]{x0394}\widehat{x}_{p}\,\unicode[STIX]{x0394}\widehat{z}_{p},\end{eqnarray}$$

where $\unicode[STIX]{x0394}\widehat{x}_{p}=5\times 10^{-3}~\text{m}$ and $\unicode[STIX]{x0394}\widehat{z}_{p}=2\times 10^{-4}~\text{m}$ are the initial distances between tracers. Because the number of tracers defining the bolus is time invariant, the bolus size $S_{b}$ is not a function of time, as is required for the incompressible flow. In figure 5(a), the dark blue cluster representing the bolus has size $5.14\times 10^{-4}~\text{m}^{2}$.

The initial shape of the bolus at $t_{0}$ is important as it describes the portion of the fluid domain transported up the slope. For the presented sample case, it is initially a thin sliver of fluid just above the boundary with an aspect ratio of approximately $60:1$. However, there is a ‘hook’ on top of the sliver, localized around $(0.05,0.003)~\text{m}$, which is also entrained into the vortex and carried upslope as presented in figure 5(b). Both the size and shape of the bolus vary with system parameters, as will be first demonstrated in § 4.1.

To contrast the properties of the boluses studied here with previous studies, we will introduce some relevant dimensionless parameters. These parameters depend on incoming wave characteristics that are computed here using the mid-depth isopycnal $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{H/2}$ at time $t_{0}$: the amplitude $a$, defined as the crest vertical displacement; the wavelength $\unicode[STIX]{x1D706}$, distance between the two zero isopycnal displacement points surrounding the crest; the wavenumber $k=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}$; and the instantaneous wave speed $c_{x}$, computed by tracking the wave crest. The typically investigated dimensionless numbers are the topographic slope $s$, the wave steepness $ka$, the internal Iribarren number $Ir=s/\sqrt{a/\unicode[STIX]{x1D706}}$, which compares the steepness of the slope and that of the incoming wave (Boegman et al. Reference Boegman, Ivey and Imberger2005; Sutherland et al. Reference Sutherland, Barrett and Ivey2013), the wave Reynolds number $Re_{w}=c_{x}ka^{2}/\unicode[STIX]{x1D708}$, the wave Richardson number $Ri_{w}=k\unicode[STIX]{x1D6FF}/(ka)^{2}$ (Thorpe Reference Thorpe1968; Troy & Koseff Reference Troy and Koseff2005) and the wave Froude number $Fr=\unicode[STIX]{x1D714}a/(g^{\prime }H/4)^{1/2}$, with $g^{\prime }=g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D70C}_{H/2}$ (Moore et al. Reference Moore, Koseff and Hult2016).

3.3 Comparison to three-dimensional analysis

While Venayagamoorthy & Fringer (Reference Venayagamoorthy and Fringer2007) demonstrated that a three-dimensional simulation yielded similar results to its two-dimensional counterpart, we must verify that the coherent structure analysis is not subject to greater coherence due to unrealistic two-dimensional turbulence (Arthur & Fringer Reference Arthur and Fringer2014). To do so, an equivalent three-dimensional simulation was performed with identical physical parameters. To perform this simulation, it was necessary to decrease the grid spatial resolution along the shelf to $\unicode[STIX]{x0394}_{x}=10^{-3}~\text{m}$. The three-dimensional domain was 0.15 m across the shelf with a spatial resolution of $\unicode[STIX]{x1D6E5}_{y}=1.5\times 10^{-3}~\text{m}$ in the $\widehat{y}$ direction, corresponding to 101 mesh layers and approximately 40 million mesh cells.

Periodic boundary conditions were applied to the lateral boundaries, and initial perturbations were incorporated to the density field in the three-dimensional simulation, to trigger lateral instabilities during the breaking. The perturbation method was inspired by Arthur & Fringer (Reference Arthur and Fringer2014, Reference Arthur and Fringer2016), and adapted to the continuously stratified system by adding the noise term $\unicode[STIX]{x1D701}^{\prime }R$, where $\unicode[STIX]{x1D701}^{\prime }=10^{-3}~\text{m}$ and $R\in [-1,1]$ is a uniformly distributed random number, to the vertical coordinate $z$ in (2.5) when initializing $\unicode[STIX]{x1D70C}_{0}(z)$, for every mesh cell. For the coherent structure analysis, the in-plane tracer spacing was not changed, and the across shelf tracer spacing was $\unicode[STIX]{x0394}\widehat{y}_{p}=0.01~\text{m}$, totalling 15 layers of tracers. For numerical purposes, the size of the grid of tracers considered in the analysis was reduced so that the total number of trajectories was $n\approx 190\,000$. The pairwise distance function between trajectories was adapted to account for the periodicity in the $\widehat{y}$ direction.

The results of the three-dimensional coherent structure analysis are presented in figure 7. A perspective of the three-dimensional bolus is presented in figure 7(a) for three time instances, $t=38$, 45 and 55 s. There is low lateral variability in the bolus shape for initial times, and only in later stages of the breaking does the development of out-of-plane motion become more pronounced (see $t=55~\text{s}$). The initial position of the bolus tracers for the three-dimensional simulation align closely to the contour representing the two-dimensional cluster as demonstrated in figure 7(b). The three-dimensional bolus is slightly bigger, with an average cross-section area of $S_{b}=5.55\times 10^{-4}~\text{m}^{2}$ as compared to the $5.14\times 10^{-4}~\text{m}^{2}$ in the two-dimensional case. In the three-dimensional case, the hook of the bolus becomes less pronounced and varies in shape for different $\widehat{x}\widehat{z}$ cross-sections. The propagation of these tracers over time allows the lateral variability to appear. The bolus tracers at $t=55~\text{s}$ are presented in figure 7(c). While the front of the three-dimensional bolus is slightly further upslope, there is little difference in the distribution of the two- and three-dimensional tracers other than a wider distribution of particle positions off the sloping boundary in the three-dimensional case. The maximum displacement upslope is approximately the same in both cases, $D_{b}=0.382~\text{m}$ for the three-dimensional case, compared to 0.402 m for the two-dimensional case.

Figure 7. (a) Time evolution of the three-dimensional bolus. The tracer lateral coordinate $\widehat{y}$ at $t_{0}=19.25~\text{s}$ is indicated in colour (from red to blue), and the different levels of shading indicate different time instances. (b) Projected two-dimensional (2-D; red) and three-dimensional (3-D; grey) bolus tracers at the initial time $t_{0}$ and (c) at $t=55~\text{s}$. Transparency is used to indicate lower or higher concentration of tracers in the projected views.

The similarities of the two- and three-dimensional results align well with the results of Venayagamoorthy & Fringer (Reference Venayagamoorthy and Fringer2007). In their case, they use a background linear stratification which results in lower reduced gravity effects, similar to the broader pycnoclines in our study. In their study, they identify lateral variation of the bolus primarily when the speed of the bolus slows to below the carrying wave speed. As the wave shoals, the water column height continues to decrease and the carrying wave and bolus decelerate at similar rates as was observed in figure 6(c). This deceleration may be the reason why there is a delayed onset of the lateral instability, which produces an alternating pattern of lobes and clefts common to the head of gravity currents flowing over a no-slip boundary (Simpson Reference Simpson1972). While there is significant across-shelf variation at successive breaking events, lateral variability is minimal for the first bolus and therefore the two-dimensional model is physically realistic, as anticipated by Aghsaee et al. (Reference Aghsaee, Boegman and Lamb2010). This limits our two-dimensional study to just the first bolus.

4.1 Pycnocline thickness variation

While past internal wave bolus studies have focused on two-layer density stratifications, the impact of the pycnocline thickness $\unicode[STIX]{x1D6FF}$ determining how smoothly the density changes with depth is investigated here. Nine different pycnocline thicknesses,

(4.1)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}=0.025,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4~\text{m},\end{eqnarray}$$

are considered. Small $\unicode[STIX]{x1D6FF}$ values represent a thin, sharp transition of density and larger values represent broad, smooth transitions. The case $\unicode[STIX]{x1D6FF}=0.025~\text{m}$ is our approximate two-layer density system. A finer transition would require a higher resolution mesh and more computational time. The other extreme pycnocline thickness case, $\unicode[STIX]{x1D6FF}=0.4~\text{m}$, matches the height of the simulation domain. The case $\unicode[STIX]{x1D6FF}=0.2~\text{m}$, presented in §§ 2 and 3, is an intermediate value between these two extremes.

While the pycnocline thickness varies, the boluses are generated by equally energetic waves, with instantaneous kinetic energy at the breaking point $E_{k,0}$ (as discussed in § 2.1). Because injecting a constant amount of energy may be simpler experimentally, that study is presented in the supplementary material. The values of $E_{k,0}$, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ and $s$ are held constant for this study. As $\unicode[STIX]{x1D6FF}$ increases, the wave amplitude increases from $a=0.012$ to 0.016 m, the wavelength remains approximately constant, $\unicode[STIX]{x1D706}=0.68\pm 0.01~\text{m}$ and the wave speed decreases from $c_{x}=0.112$ to $0.088~\text{m}~\text{s}^{-1}$. The wave steepness $ka$ varies from 0.11 to 0.15, the internal Iribarren number $Ir$ from 1.3 to 1.1, the wave Reynolds number $Re_{w}$ from 160 to 225, the Richardson number $Ri_{w}$ from 17 to 157 and the Froude number $Fr$ from 0.056 to 0.074, all of them monotonically with $\unicode[STIX]{x1D6FF}$.

Figure 8. Time evolution of the boluses from $t_{0}$ to $t_{f}$ for stratifications with pycnocline thicknesses (a) $\unicode[STIX]{x1D6FF}=0.025~\text{m}$, (b) 0.1 m, (c) 0.2 m, (d) 0.25 m and (e) 0.3 m. The trajectory of the bolus centre of volume for $\unicode[STIX]{x1D6FF}=0.2~\text{m}$ and the displacement upslope $D_{b}$ are illustrated in (c). (Video of the evolution of the figure is available in the supplementary material.)

The impact of the pycnocline thickness on the bolus dynamics is presented in figure 8. Snapshots of the bolus in the rotated frame $\widehat{O}_{\widehat{x}\widehat{z}}$ are presented for selected times within the respective time span $[t_{0},t_{f}]$. For the sample case, in figure 8(c), the trajectory of the centre of volume is presented in grey and the maximum displacement upslope $D_{b}$ defined in (3.4) is illustrated. The values of $D_{b}$ and the bolus lifetime ($t_{f}-t_{0}$) both grow monotonically with $\unicode[STIX]{x1D6FF}$. A more complex relationship is observed between the bolus size $S_{b}$ and $\unicode[STIX]{x1D6FF}$: boluses are smaller for the approximate two-layer system, their sizes initially increase with $\unicode[STIX]{x1D6FF}$, reach a maximum, then shrink from this maximum as $\unicode[STIX]{x1D6FF}$ approaches the broadest pycnocline case. It is also worth noting that as $\unicode[STIX]{x1D6FF}$ increases the wave breaking happens later in time and the fluid trapped is initially located further up the slope.

Figure 9. (a) Bolus size $S_{b}$ and (b) maximum displacement upslope $D_{b}$ as a function of the pycnocline thickness, $\unicode[STIX]{x1D6FF}$, for mode-1 waves breaking with constant energy. The error bars represent the range of displacement within the bolus tracers. Different markers are used to represent bolus shape categories: ball (circle), hook (triangle) and sliver (square).

The bolus size $S_{b}$ and displacement upslope $D_{b}$ are presented in figure 9. The maximum $S_{b}$ is observed for $\unicode[STIX]{x1D6FF}=0.15~\text{m}$, while $D_{b}$ monotonically increases with $\unicode[STIX]{x1D6FF}$. The largest bolus is approximately twice the size of the one observed for the approximate two-layer simulation and five times as large as the bolus produced in the broadest pycnocline stratification. The bolus obtained for the broadest pycnocline case, however small, travels about three times as far as the two-layer bolus, with a lifetime more than twice as long. Upper and lower bounds for the range of tracer displacement within the bolus are also presented in figure 9(b). This range is smallest for the narrowest pycnoclines and gradually grows until $\unicode[STIX]{x1D6FF}=0.3~\text{m}$. Both the upper and lower bound increase monotonically.

Figure 10. (a) Ball, (b) hook and (c) sliver bolus shape categories represented by the bolus tracer position at time $t_{0}$. Shapes are obtained for $\unicode[STIX]{x1D6FF}=0.025~\text{m}$ (blue), 0.25 m (green) and 0.3 m (red). The aspect ratio is approximately 30:1. The markers in the top right of each plot are used to represent bolus shape categories.

A jump in $S_{b}$ and $D_{b}$ occurs between boluses formed with $\unicode[STIX]{x1D6FF}=0.25~\text{m}$ and $\unicode[STIX]{x1D6FF}=0.3~\text{m}$. This result relates to a difference in the shape of the bolus for each case. We refer to the bolus shape based on the initial geometry of the bolus tracers. In this paper, bolus shapes are divided into three categories: ball, hook and sliver, which are illustrated in figure 10. The ball category, illustrated in figure 10(a), is characterized by having a convex geometry, and is typical of thin pycnoclines, for which the breaking is more abrupt. Increasing $\unicode[STIX]{x1D6FF}$ causes the initial bolus shape to take on the shape of a ‘hook’, which is eventually entrained in the vortex during the bolus propagation. The hook category is well illustrated by the case $\unicode[STIX]{x1D6FF}=0.25~\text{m}$ in figure 10(b), in which the hook is located around $(0.1,0.003)~\text{m}$. As $\unicode[STIX]{x1D6FF}$ grows further, only a horizontal stripe of tracers remain part of the bolus, and the hook is lost. The single horizontal stripe pattern defines the sliver category, illustrated in figure 10(c). The loss of the hook is a result of the propagating vortex not being strong enough to trap the material located in that area. In fact, the maximum vorticity felt by the bolus monotonically decreases as $\unicode[STIX]{x1D6FF}$ increases. The transition from hook to sliver occurs somewhere between $\unicode[STIX]{x1D6FF}=0.25$ and $\unicode[STIX]{x1D6FF}=0.3~\text{m}$. This transition results in the bolus size decreasing by more than half. This significant decrease in size may be the reason for the smaller range in displacements observed for $\unicode[STIX]{x1D6FF}=0.3~\text{m}$ compared to $\unicode[STIX]{x1D6FF}=0.25~\text{m}$. The corresponding increase in $D_{b}$ is related to the fact that the tracers that form the hook component are typically the first ones to be ejected from the vortex.

Drastic changes in the bolus size and displacement are found to be related to transitions in the bolus geometric properties described by these three categories. For this reason, for figure 9 and all the $S_{b}$ and $D_{b}$ plots in the following sections, distinct bolus shapes are plotted using different symbols (as introduced in figure 10), so that transitions in shape are easily identifiable.

4.2 Shoaling wave energy variation

In the previous section, waves reached the slope with equal kinetic energy $E_{k,0}$. Here, the influence of the shoaling wave energy on the bolus properties is investigated. We adjust the velocity forcing amplitude $a_{f}$ to produce waves with the energy scale factors

(4.2)$$\begin{eqnarray}E_{k}/E_{k,0}=1/8,1/4,1/2,1,2,4.\end{eqnarray}$$

The upper bound for energy scale factors considered is set by computational limitations of direct simulations in highly energetic cases. The lower bound attempts to capture low-energy cases in which boluses do not form. To understand how the wave energy impacts the bolus propagation as a function of the pycnocline thickness, three thicknesses ($\unicode[STIX]{x1D6FF}=0.025,0.2$ and $0.4~\text{m}$) are considered for each energy scale factor. For all simulations, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ and $s$ are held constant. The wave amplitude increases with energy from $a=0.004$ to 0.033 m, $\unicode[STIX]{x1D706}$ remains approximately constant and $c_{x}$ is minimally impacted by the energy scale factor, being mostly determined by $\unicode[STIX]{x1D6FF}$. Highly energetic waves are therefore steeper ($ka=0.04$ to 0.30), resulting in $Ir$ varying from 2.19 to 0.80. The value of $Re_{w}$ grows from 20 to 855, $Ri_{w}$ decreases with $E_{k}$ and grows with $\unicode[STIX]{x1D6FF}$, spanning 6 to 1200, and $Fr$ grows with both $E_{k}$ and $\unicode[STIX]{x1D6FF}$ from 0.02 to 0.15.

Figure 11. (a) Bolus size $S_{b}$ and (b) maximum displacement upslope $D_{b}$ as a function of the energy factor for the shoaling wave. Simulations were performed with the pycnocline thicknesses $\unicode[STIX]{x1D6FF}=0.025$ (blue), 0.2 (black) and 0.4 m (red). When no bolus is identified, a cross is plotted at $S_{b}=0~\text{m}^{2}$ and $D_{b}=0~\text{m}$. Marker shapes are used to represent each bolus shape categories: ball (circle), hook (triangle) and sliver (square).

Figure 11 presents the relationship between the shoaling wave energy and the bolus properties. The bolus size dependence is presented in figure 11(a) and the displacement dependence in figure 11(b). When no bolus is identified, the properties for that case are plotted in the graph as crosses with $S_{b}=0~\text{m}^{2}$, $D_{b}=0~\text{m}$. When the energy scale factor is $1/8$, none of the three pycnocline thicknesses produce a bolus. For $E_{k}/E_{k,0}=1/4$, the only case producing a bolus is $\unicode[STIX]{x1D6FF}=0.2~\text{m}$. For $E_{k}/E_{k,0}\geqslant 1/2$, all three cases of $\unicode[STIX]{x1D6FF}$ produce boluses, and $S_{b}$ increases with the wave energy. The bolus size sensitivity to $E_{k}$ is weakest for the two-layer case, which is possibly related to the more abrupt breaking and quick mixing characteristic of this stratification. The pycnocline thickness producing the largest and smallest bolus depends on $E_{k}$.

Varying the wave energy also results in transitioning between bolus shape categories, and transitions happen between different energy factors for different pycnocline thicknesses. In particular, for $\unicode[STIX]{x1D6FF}=0.2~\text{m}$ a transition from sliver to hook-shaped boluses occurs when the energy factor is increased from $1/2$ to 1, while this same transition for $\unicode[STIX]{x1D6FF}=0.4~\text{m}$ is observed between the energy factors 1 to 2. This result is related to the fact that, at higher energy, the vortex has a greater circulation and is more capable of retaining particles that are moving with the bolus. The jumps in $S_{b}$ and $D_{b}$ resulting from these shape transitions produce counter-intuitive results: for example, the bolus displacement for $\unicode[STIX]{x1D6FF}=0.2~\text{m}$ is larger when the incoming wave has half of the energy compared to the baseline energy. The energy factor $1/2$ produces a sliver bolus that is approximately half the size of the hook bolus of energy $E_{k,0}$, so this bolus consists of fewer tracers that travel a longer distance. Such hook–sliver transitions explain why there are two clearly different growth rates observed for $S_{b}$ and $D_{b}$ for $\unicode[STIX]{x1D6FF}=0.2$ and 0.4 m, as the smaller sliver-shaped boluses travel longer distances compared to the larger hook-shaped boluses.

4.3 Density change variation

The stratification was varied in the previous sections by changing the pycnocline thickness, while keeping the density change parameter $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ constant. Variations in $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ have not been investigated in previous bolus studies. However, when considering typical ocean stratification profiles, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ varies and its value may have an impact on the resulting bolus dynamics. Here, we vary the stratification by imposing density changes

(4.3)$$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=10,20,40,80~\text{kg}~\text{m}^{-3}.\end{eqnarray}$$

The choice of not studying values greater than $80~\text{kg}~\text{m}^{-3}$ is related to the limitations of the Boussinesq approximation when describing stronger density variations. For the lower bound, smaller changes than $10~\text{kg}~\text{m}^{-3}$, for the fixed value of $\unicode[STIX]{x1D714}$, would result in a stratification with buoyancy frequency $N(z)<\unicode[STIX]{x1D714}$ for all $z$, resulting in evanescent waves. To understand how the sensitivity to $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ varies across different pycnocline thicknesses, each density change is studied for the values of $\unicode[STIX]{x1D6FF}$ from 0.025 to 0.25 m. In all cases, the energy of the shoaling wave is $E_{k,0}$ and the topographic slope $s$ is held constant.

The combination of different density changes and pycnocline thicknesses has a significant impact on all of the wave characteristic: $a$ varies from 0.023 to 0.004 m, $\unicode[STIX]{x1D706}$ from 0.65 to 1.21 m and $c_{x}$ from 0.08 to $0.25~\text{m}~\text{s}^{-1}$. For higher values of $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$, the waves have smaller amplitude and larger wavelength and speed. For this study, the wave steepness varies within $ka=0.04$ and 0.22, so that $Ir=0.95$ to 3.06. The value of $Re_{w}$ decreases with $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$, from 378 to 22, $Ri_{w}$ increases from 7 to 1980 while $Fr$ decreases from 0.15 to 0.01.

Figure 12. (a) Bolus size $S_{b}$ and (b) maximum displacement upslope $D_{b}$ as a function of the pycnocline thickness, for variable density change, $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$. Simulations were performed with density changes $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=10$ (red), 20 (black), 40 (blue) and $80~\text{kg}~\text{m}^{-3}$ (green). When no bolus is identified, a cross is plotted at $S_{b}=0~\text{m}^{2}$ and $D_{b}=0~\text{m}$. Marker shapes are used to represent each bolus shape categories: ball (circle), hook (triangle) and sliver (square).

The bolus size and displacement upslope are presented in figure 12. Figure 12(a) demonstrates that the pycnocline thickness that maximizes the bolus size depends on $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$. The pycnocline thickness that results in the largest bolus increases from $\unicode[STIX]{x1D6FF}=0.1~\text{m}$ for $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=80~\text{kg}~\text{m}^{-3}$ to $\unicode[STIX]{x1D6FF}=0.2~\text{m}$ for $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=10~\text{kg}~\text{m}^{-3}$. For a strongly stratified system, the wave essentially reflects from the sloping wall, forming either a small bolus or no bolus, as is the case for $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=80~\text{kg}~\text{m}^{-3}$ and $\unicode[STIX]{x1D6FF}=0.025~\text{m}$. This case corresponds to one the greatest internal Iribarren numbers observed, $Ir=2.97$. Figure 12(b) confirms the trend previously observed in figure 9(b): the displacement $D_{b}$ increases monotonically with $\unicode[STIX]{x1D6FF}$ for every $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ value considered.

One of the factors impacting $S_{b}$ for different stratifications is the speed of the internal wave, which is comparable to the speed of the bolus upslope (see figure 6). Based on our simulations, slowly propagating boluses tend to trap more material inside the vortex than faster boluses. Another factor that plays a role in how much material is transported is the dynamics of the breaking itself. Larger buoyancy frequencies result in the bolus edge destabilizing and mixing more quickly, resulting in less effective transport.

4.4 Topographic slope variation

Finally, to investigate the influence of the underlying topography, the topographic slope $s$, held at 0.176 for the previous studies, is varied. Keeping a constant density change $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}=20~\text{kg}~\text{m}^{-3}$ and shoaling wave energy $E_{k,0}$, the topographic slopes considered are

(4.4)$$\begin{eqnarray}s=0.105,0.123,0.141,0.158,0.176,0.194,0.213,0.231,\end{eqnarray}$$

which correspond to slope angles of $6^{\circ }$–$13^{\circ }$. To ensure that the wave travels the same distance $L$ before breaking on slope, the domain geometry is adjusted so that the horizontal plane $z=H/2$ intersects the topography at $x=L$, as presented in figure 2. Such a geometric constraint increases the domain size (and thus the number of mesh cells) for smaller slopes. The tested slopes are also constrained so that the topography does not intersect the domain inlet. For these two reasons $s$ is not decreased below 0.105. To identify the relationship between bolus transport, topographic slope and pycnocline thickness, the thicknesses $\unicode[STIX]{x1D6FF}=0.025,0.2$ and 0.4 m are tested with each of the slopes. Three slopes corresponding to the extremes and the baseline case ($s=0.105,0.176$ and 0.231) are also tested with all nine pycnocline thicknesses. While varying the topographic slope $s$ does not alter the properties of the incoming wave nor the dimensionless parameters $Re_{w}$, $Ri_{w}$ or $Fr$ from § 4.1, changes in $s$ affect the internal Iribarren number $Ir$, which spans the range 0.67 through 1.71 as $s$ is increased.

Figure 13. (a) Bolus size $S_{b}$, (b) maximum displacement upslope $D_{b}$ and (c) maximum vertical displacement $\unicode[STIX]{x0394}z_{b}$ as a function of the topographic slope $s$ for the pycnocline thicknesses $\unicode[STIX]{x1D6FF}=0.025$ (blue), 0.2 (black), 0.4 m (red). (d) Bolus size $S_{b}$, (e) maximum displacement upslope $D_{b}$ and (f) maximum vertical displacement $\unicode[STIX]{x0394}z_{b}$ as a function of the pycnocline thickness $\unicode[STIX]{x1D6FF}$, for slopes $s=0.105$ (red), 0.176 (black), 0.231 (blue). Marker shapes are used to represent each bolus shape category: ball (circle), hook (triangle) and sliver (square). Empty red markers for $s=0.105$, $\unicode[STIX]{x1D6FF}=0.3~\text{m}$ represent two possible values for the bolus at the hook-sliver transition.

The slope variation results are presented in figure 13. Initially, we examine the direct dependence on the slope $s$ for the three values of $\unicode[STIX]{x1D6FF}$ (figure 13a–c). The dependence of the bolus size $S_{b}$ on $s$ is presented in figure 13(a). For $\unicode[STIX]{x1D6FF}=0.025$ and 0.4 m, maxima in bolus size are obtained at approximately $s=0.123$. For $\unicode[STIX]{x1D6FF}=0.2~\text{m}$, there is a monotonic decrease in the bolus size as the topographic slope is increased. Based on the behaviour observed for the extreme cases, there must be a slope shallower than $s=0.105$ that maximizes the size of the bolus for the case $\unicode[STIX]{x1D6FF}=0.2~\text{m}$ because there is ultimately no bolus generated when there is no slope. With respect to the displacement upslope $D_{b}$, presented in figure 13(b), steeper slopes result in boluses whose displacement along the slope decreases monotonically with $s$. This same trend is not observed if, instead of the displacement upslope $D_{b}$, the quantity plotted is the effective vertical displacement $\unicode[STIX]{x0394}z_{b}$. While $D_{b}$ decreases with $s$, $\unicode[STIX]{x0394}z_{b}$ is observed to remain nearly constant as $s$ is varied. This nearly constant $\unicode[STIX]{x0394}z_{b}$ behaviour indicates that the stratification alone plays a major role in setting $\unicode[STIX]{x0394}z_{b}$, and $D_{b}$ is the projection of $\unicode[STIX]{x0394}z_{b}$ along the slope, resulting in longer displacements for shallower slopes and weaker stratifications.

Finally, the direct dependence of $S_{b}$, $D_{b}$ and $\unicode[STIX]{x0394}z_{b}$ on $\unicode[STIX]{x1D6FF}$ for three values of the slope $s$ is presented in figure 13(d–f). The bolus size dependence in figure 13(d) demonstrates that, for all three values of $s$, the maximum $S_{b}$ is obtained for $\unicode[STIX]{x1D6FF}=0.15~\text{m}$. That means that the stratification corresponding to maximum bolus size is independent of the topographic slope for the range of values considered. The displacement upslope is presented in figure 13(e) and reveals that the displacement increases nearly monotonically with $\unicode[STIX]{x1D6FF}$. Figure 13(f) presents the corresponding vertical displacement $\unicode[STIX]{x0394}z_{b}$, and a close correspondence between the different $s$ cases is observed.

It is important to notice that the topographic slope also plays a role in determining when the ball–hook or hook–sliver bolus shape category transitions occur with respect to the pycnocline thicknesses. Also, note that the combination $s=0.105$ and $\unicode[STIX]{x1D6FF}=0.3~\text{m}$ leads to the single instance in this paper where the bolus is very close to the transitioning point between hook and sliver categories, and both outputs are plotted in figure 13(d–f) with empty symbols. The output from the method as described in § 3.1 corresponds to the hook case, but when we verify the robustness of the analysis by adding an extra cluster in the K-means step, the analysis produces a sliver-shaped bolus. This was the only case where this robustness analysis modified the bolus.

4.5 Dimensionless analysis

The parametric studies presented in §§ 4.1–4.4 focus on varying one physical parameter at a time and evaluating how that parameter, in conjunction with the pycnocline thickness, impacts the properties of the bolus. Previous efforts have utilized dimensionless parameters to provide arguments about how the bolus will break, properties of the bolus and the amount of mixing induced by the breaking process (Venayagamoorthy & Fringer Reference Venayagamoorthy and Fringer2007; Aghsaee et al. Reference Aghsaee, Boegman and Lamb2010; Sutherland et al. Reference Sutherland, Barrett and Ivey2013; Moore et al. Reference Moore, Koseff and Hult2016; Arthur et al. Reference Arthur, Koseff and Fringer2017). Most of these studies focused on shoaling internal solitary waves of depression in a two-layer or nearly two-layer density stratification, where the forcing mechanism and the resulting wave properties can be directly controlled. Because the mode-1 generated, nonlinear internal waves simulated here experience some dispersion as they propagate from the inlet, studying the isolated effects of the dimensionless numbers or the wave properties is particularly challenging. Nevertheless, it is possible to recast our physical studies in a dimensionless perspective to relate our findings to previous studies.

Figure 14. Relationship between dimensionless parameters and bolus size $S_{b}$ and distance upslope $D_{b}$. Marker sizes represent the magnitude of the plotted quantities, and colours represent different studies: pycnocline thickness variation (grey), energy variation (green), density change variation (blue), slope variation (red). Crosses represent cases for which no bolus is formed.

The size and displacement dependence on pairs of dimensionless values defined in § 3.2 is presented in figure 14. The first comparison evaluates the relationship between $Ir$ and $k\unicode[STIX]{x1D6FF}$ in figure 14(a,b). Here, the size of the bolus decreases with increasing $Ir$ and there is an optimal $k\unicode[STIX]{x1D6FF}$ that produces the largest bolus. The optimal $k\unicode[STIX]{x1D6FF}$ appears to decrease as $Ir$ increases. The displacement of the bolus increases for increasing $k\unicode[STIX]{x1D6FF}$ and decreasing $Ir$. This trend is observed in figure 14(c,d) as $a$ is correlated to $\unicode[STIX]{x1D6FF}$. This results in an optimal $ka$ for bolus size. The displacement upslope increases with $ka$ and decreases as the topographic slope increases. However, the trend in figure 13(c) indicates that the vertical displacement of the bolus has a more complicated relationship with the slope. Next, $Re_{w}$ and $Ri_{w}$ are compared in figure 14(e,f). These two numbers are most strongly varied by changing the energy at the breaking location $E_{k}$. Based on the simulations performed, the size of the bolus increases with $Re_{w}$, but there is a non-monotonic trend with $Ri_{w}$. The bolus displacement grows with both $Re_{w}$ and $Ri_{w}$. It is interesting to note that for $Re_{w}<80$ no bolus is produced. Finally, the dependence with $Fr$ is presented in figure 14(g,h). These figures demonstrate that the bolus size grows with $Fr$ but there is limited change in displacement.

The calculation of these dimensionless numbers allows us to relate the bolus breaking to those observed in other studies. Aghsaee et al. (Reference Aghsaee, Boegman and Lamb2010) numerically studied the breaking of two-layer waves for a range of topographic and wave slopes. Our topographic slope ranges from $s=0.105$ to 0.231 and the wave slope $a/\unicode[STIX]{x1D706}$ ranges from 0.017 to 0.025. Based on their study, the corresponding wave breakers will be classified as ‘surging’ events, which are characterized by minimal breaking of the incoming wave interface, surging of the wave upslope and the production of a bolus. This aligns well with what is observed in our simulations. Arthur & Fringer (Reference Arthur and Fringer2016) also simulated the breaking of two-layer waves, observing surging breakers for their waves of lowest amplitude (case 1), which still have higher amplitude than any of the waves addressed here. Finally, Arthur et al. (Reference Arthur, Koseff and Fringer2017) considered higher-amplitude waves for different pycnocline thicknesses and classified all of their breakers as ‘collapsing’, which falls within the appropriate two-layer breaker classification, despite the pycnocline thickness variation. This matches our observation that changing the pycnocline thickness does not seem to modify the breaking classification.