1. Introduction
Internal solitary waves (ISWs) are a fundamental phenomenon in stably stratified fluids such as oceans, lakes, estuaries and the atmosphere. These waves are driven by gravitational restoring forces acting on vertically displaced fluid parcels and propagate horizontally along density interfaces, typically the pycnocline (where density varies quickly with depth). The ISWs are nonlinear and dispersive and can travel long distances without a significant change in shape. In the ocean, their amplitudes, as measured by the maximum wave-induced isopycnal displacement, are a significant fraction of the water column depth (Duda et al. Reference Duda, Lynch, Irish, Beardsley, Ramp, Chiu, Tang and Yang2004; Huang et al. Reference Huang, Chen, Zhao, Zhang, Zhou, Yang and Tian2016), making them a key mechanism for transferring energy from large-scale forcing (e.g. tides and wind) to small-scale turbulent dissipation and mixing. The ISWs have wide-ranging implications for ocean circulation, climate modelling, marine biogeochemistry, sediment resuspension and transport, underwater communication and offshore engineering. Due to their complex dynamics and environmental significance, ISWs have been the focus of extensive research. Excellent reviews are available on: (i) nonlinear internal waves (Helfrich & Melville Reference Helfrich and Melville2006), (ii) mechanisms of ISW generation (Jackson et al. Reference Jackson, Da Silva and Jeans2012), (iii) breaking and dissipation of ISW on continental slopes (Lamb Reference Lamb2014; Lamb et al. Reference Lamb, Lien and Diamessis2019), (iv) resuspension and transport of sediments due to ISWs (Boegman & Stastna Reference Boegman and Stastna2019) and (v) shoaling of ISWs (Stastna & Legare Reference Stastna and Legare2024).
As ISWs shoal (e.g. on continental slopes and lake shores), they undergo significant transformation as a result of increasing interaction with the bottom boundary as the water depth decreases. This includes wave breaking and energy dissipation in both the bottom boundary layer (BBL) and the pycnocline. These processes are central to the global energy cascade from large-scale mechanical forcing to small-scale mixing (St. Laurent et al. Reference St. Laurent, Simmons, Tang and Wang2011; Sarkar & Scotti Reference Sarkar and Scotti2017). Shoaling ISWs can also upwell cold, nutrient-rich water, which benefits benthic ecosystems such as coral reefs by improving their resilience to climate change (Green et al. Reference Green, Jones, Rayson, Lowe, Bluteau and Ivey2019; Reid et al. Reference Reid, DeCarlo, Cohen, Wong, Lentz, Safaie, Hall and Davis2019; Wyatt et al. Reference Wyatt, Leichter, Toth, Miyajima, Aronson and Nagata2020; Pérez-Rosales et al. Reference Pérez-Rosales, Rouzé, Torda, Bongaerts, Pichon, Parravicini and Hédouin2021). In addition, shoaling is associated with: (i) instabilities and turbulence in the BBL (Bogucki et al. Reference Bogucki, Dickey and Redekopp1997; Chang et al. Reference Chang2021; Trowbridge et al. Reference Trowbridge, Helfrich, Reeder, Medley, Chang, Jan, Ramp and Yang2025), (ii) enhanced sediment resuspension and transport (Hosegood et al. Reference Hosegood, Bonnin and van Haren2004; Quaresma et al. Reference Quaresma, Vitorino, Oliveira and da Silva2007; Bourgault et al. Reference Bourgault, Morsilli, Richards, Neumeier and Kelley2014; Zulberti et al. Reference Zulberti, Jones and Ivey2020; Trowbridge et al. Reference Trowbridge, Helfrich, Reeder, Medley, Chang, Jan, Ramp and Yang2025), (iii) nutrient cycling in the euphotic zone (Sandstrom & Elliott Reference Sandstrom and Elliott1984) and (iv) the formation of sedimentary structures such as subaqueous sand dunes (Reeder et al. Reference Reeder, Ma and Yang2011; Droghei et al. Reference Droghei, Falcini, Casalbore, Martorelli, Mosetti, Sannino, Santoleri and Chiocci2016; Ghassemi & Boegman Reference Ghassemi and Boegman2024).
The most commonly observed ISWs in nature are mode-1, which displace isopycnals in a single direction, either downward (depression) or upward (elevation). The wave form depends on the pycnocline’s position relative to mid-depth: depression waves occur when the pycnocline lies above mid-depth, while elevation waves occur when it lies below. As mode-1 depression ISWs shoal on gentle slopes, they fission into a train of elevation waves, as observed by Orr & Mignerey (Reference Orr and Mignerey2003), Zhao et al. (Reference Zhao, Klemas, Zheng, Li and Yan2004), Shroyer et al. (Reference Shroyer, Moum and Nash2009), St. Laurent et al. (Reference St. Laurent, Simmons, Tang and Wang2011) and Sinnett et al. (Reference Sinnett, Ramp, Yang, Chang, Jan and Davis2022). Recent observations (Sinnett et al. Reference Sinnett, Ramp, Yang, Chang, Jan and Davis2022; Whitwell et al. Reference Whitwell, Jones, Ivey, Rosevear and Rayson2024) also highlight the role of elevation ISWs, mode-2 ISWs and internal boluses in shelf sea mixing and sediment resuspension during shoaling events.
Internal solitary wave shoaling can occur with and without slope-intersecting pycnoclines. This study will focus on the former case, for the latter case, interested readers are referred to the two-dimensional simulations of Lamb & Warn-Varnas (Reference Lamb and Warn-Varnas2015) motivated by the South China Sea, field studies of Lien et al. (Reference Lien, D’Asaro, Henyey, Chang, Tang and Yang2012, Reference Lien, Henyey, Ma and Yang2014) focussing on a wave with a sub-surface trapped core and the recent work of Bolioudakis et al. (Reference Bolioudakis, Diamantopoulos, Diamessis, Lien, Lamb, Rivera-Rosario and Thomsen2026) which represents three-dimensional simulations for this type of wave, revealing convective and shear instabilities within the core of the ISW as it shoals.
Numerous modelling studies have examined the shoaling of mode-1 depression ISWs over smooth, idealised, linear topographic slopes, where the pycnocline intersects the slope, using laboratory experiments (e.g. Michallet & Ivey Reference Michallet and Ivey1999; Boegman et al. Reference Boegman, Ivey and Imberger2005; Sutherland et al. Reference Sutherland, Barrett and Ivey2013) and numerical simulations (e.g. Aghsaee et al. Reference Aghsaee, Boegman and Lamb2010; Arthur & Fringer Reference Arthur and Fringer2014; Nakayama et al. Reference Nakayama, Sato, Shimizu and Boegman2019). These works have established a classification of breaking types, namely plunging, collapsing, surging and fission, based on wave and slope steepness. Analogous to surface wave studies, an internal Iribarren number (wave steepness to slope steepness ratio) is used to categorise breaking behaviour (Boegman et al. Reference Boegman, Ivey and Imberger2005; Sutherland et al. Reference Sutherland, Barrett and Ivey2013). Nakayama et al. (Reference Nakayama, Sato, Shimizu and Boegman2019) have further shown that an internal Reynolds number influences the classification and Hartharn-Evans et al. (Reference Hartharn-Evans, Carr, Stastna and Davies2022) demonstrated that the background stratification profile significantly influences the classification of the shoaling dynamics. Recent studies have explored particle transport and sediment resuspension (Deepwell et al. Reference Deepwell, Sapède, Buchart, Swaters and Sutherland2020), cross-boundary-layer transport (Xu & Stastna Reference Xu and Stastna2020), mixing efficiency (la Forgia et al. Reference la Forgia, Cavaliere, Adduce and Falcini2021), dense pulses from fissioning (Hartharn-Evans et al. Reference Hartharn-Evans, Stastna and Carr2023), parameterisations for large-scale models (He et al. Reference He, Chen, Li, Xu and Meng2023) and three-dimensional (3-D) simulations of BBL instabilities (Castro-Folker & Stastna Reference Castro-Folker and Stastna2024).
Given the complexity of the problem, most modelling studies have focused on smooth, linear, topography to understand the underlying bulk flow features. However, natural slopes are rarely linear or smooth and can exhibit a wide range of roughness, from sand ripples to coral reefs, rocks and kelp forests, each with distinct length scales and hydrodynamic impacts. Few studies have addressed the influence of variable topography on ISW shoaling, despite field observations (Reeder et al. Reference Reeder, Ma and Yang2011; Ramp et al. Reference Ramp, Yang, Chiu, Reeder and Bahr2021) and field-scale simulations (Lamb & Warn-Varnas Reference Lamb and Warn-Varnas2015; Rivera-Rosario et al. Reference Rivera-Rosario, Diamessis, Lien, Lamb and Thomsen2020, Reference Rivera-Rosario, Diamessis, Lien, Lamb and Thomsen2023; Zang et al. Reference Zang, Zhang, Chen, Xie, Zou and Li2023) revealing a complex dynamic over undulating slopes. Ghassemi & Boegman (Reference Ghassemi and Boegman2024) showed in laboratory experiments that sand ripples can form due to shoaling ISWs over a movable sedimentary bed. Motivated by coral reef environments, Yuan et al. (Reference Yuan, Luo, Wang, Lin, Liu and Song2024) conducted laboratory experiments on rough slopes and found that increased roughness enhances turbulent dissipation but can inhibit mixing efficiency and reduce wave energy loss. However, their study was limited to a steep slope (
$s = 0.577$
, where
$s$
is the vertical displacement per unit horizontal displacement), which restricted their shoaling wave dynamics to collapsing–surging type only. Wang et al. (Reference Wang, Yuan, Lin and Song2025) have recently investigated the effect of ISW breaking on reef-lagoon systems using laboratory experiments. They used LEGO blocks of varying height as a proxy for coral reefs of varying roughness. Like Yuan et al. (Reference Yuan, Luo, Wang, Lin, Liu and Song2024) their study was limited to ISWs of collapsing–surging type only. Stastna & Legare (Reference Stastna and Legare2024) presented illustrative numerical examples of wave shoaling over undulating topography but did not provide a systematic investigation and Posada-Bedoya et al. (Reference Posada-Bedoya, Olsthoorn and Boegman2025) have recently considered the effect of corrugated topography on BBL instabilities under ISWs but only for a horizontal bed.
This study presents, for the first time, a comprehensive numerical investigation of ISW shoaling over linear slopes with undulating topography, under the unifying condition of a slope-intersecting pycnocline, validated against previous experiments on a smooth slope (Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022) and over a horizontal corrugated bed (Carr et al. Reference Carr, Stastna and Davies2010). All four breaking types (collapsing, plunging, surging and fission) are examined. The slope steepness range considered is
$0.033 \lt s \lt 0.2$
, which includes typical values for continental slopes (
$0.03 \lessapprox s \lessapprox 0.07$
) and lakes (
$0.01 \lessapprox s \lessapprox 0.1$
). Wave steepness is varied accordingly to capture the full spectrum of the shoaling dynamics. The amplitude and frequency of the corrugations on the slope are varied. In all but the fission case the horizontal length scales of the corrugations are significantly less than the wavelength. In the fission case, the topographic slope is shallower and longer, hence a broader range of topographic wavelengths can be considered including some that exceed the ISW wavelength. To minimise the variables under consideration, the background stratification and Reynolds number are held fixed and attention is focussed on 2-D simulation. A limited number of 3-D simulations are presented to show that the onset of breaking is accurately captured by 2-D simulations and in order to aid discussion on and motivate future study around turbulence and mixing. The paper is organised as follows: § 2 details the numerical model; § 3 presents the results from the 2-D simulations; § 4 discusses 3-D simulations and what can be inferred from them; § 5 presents the user-defined scatter plot (USP) method; § 6 offers conclusions and directions for future work.
2. Numerical method
2.1. Governing equations and numerical domain
A spectral parallel incompressible Navier–Stokes code (Subich et al. Reference Subich, Lamb and Stastna2013) was utilised to solve the stratified Navier–Stokes equations subject to the Boussinesq approximation
\begin{align} {{\partial {\bf{u}}} \over {\partial t}} + {\bf{u}} \cdot \nabla {\bf{u}} &= - {1 \over {{\rho _0}}}\nabla \! P + \nu {\nabla ^2}{\bf{u}} - {{\rho g} \over {{\rho _0}}}\hat k, \\[-12pt]\nonumber \end{align}
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} & = 0, \\[-12pt]\nonumber \end{align}
\begin{align} {{\partial \rho } \over {\partial t}} + {\bf{u}} \cdot \nabla \rho &= \kappa {\nabla ^2}\rho, \end{align}
where
$\boldsymbol{u}$
is the velocity,
$P$
is the pressure,
$\rho$
is the density and
$\rho _0$
is some reference density of the fluid (here taken to be the density in the lower layer of the stratified water column). The physical parameters are acceleration due to gravity
$g=9.81\,\textrm {ms}^{-2}$
, shear viscosity
$\nu =10^{-6}\,\textrm {m}^2\,\textrm {s}^{-1}$
and scalar diffusivity
$\kappa =10^{-7}\,\textrm {m}^2\,\textrm {s}^{-1}$
. The unit vector in the vertical direction is denoted by
$\hat {k}$
. The Reynolds numbers associated with the flow (after wave generation and prior to interaction with the slope) were
$\textit{Re}_{\textit{ISW}}=O(10^2)$
,
$\textit{Re}_{w0}=O(10^4)$
and
$\textit{Re}_{w}=O(10^6)$
where
\begin{align} {\textit{Re}_{\textit{ISW}}} = {U_2}\sqrt {{\lambda \over {\nu ({U_2} + c)}}}, \quad {\textit{Re}_{w0}} = {{{c_0}H} \over \nu },\quad {\textit{Re}_w} = {{cH} \over \nu }, \end{align}
and
$U_2$
is the absolute value of the maximum horizontal velocity at the wave trough,
$\lambda$
is the wavelength at half-amplitude,
$c$
is wave speed,
$c_0$
is the linear long wave speed and
$H$
is the total fluid depth. Here,
$\textit{Re}_{\textit{ISW}}$
is an estimate of the momentum thickness Reynolds number introduced by Aghsaee et al. (Reference Aghsaee, Boegman, Diamessis and Lamb2012).
Numerical domain and typical initial stratification. Colour scheme depicts the density anomaly
$\rho '=\rho -\rho _1$
, where
$\rho _1$
is the density in the top layer and
$\rho$
is density. The topographic slope is coloured white.
The governing equations were solved on a numerical domain like that depicted in figure 1. The origin is chosen to be aligned with the toe of the topographic slope in the horizontal
$x$
-direction and at the top of the water column in the vertical
$z$
-direction. The stratification in the main part of the domain was modelled using a vertical density profile given by
\begin{equation} \rho (z) = {\rho _2} + \Delta \rho \tanh \left( {{{z - {z_{\textit{pyc}}}} \over {{h_{\textit{pyc}}}}}} \right)\!, \end{equation}
where
$\rho _2$
is the density of the lower layer,
$\Delta \rho = \rho _2-\rho _1$
is the change in density through the water column,
$\rho _1$
is the density in the top layer,
$z_{\textit{pyc}} = -0.07\,\textrm {m}$
is the depth of the centre of the pycnocline and
$h_{\textit{pyc}} = 0.003\,\textrm {m}$
is the pycnocline half-thickness. This background stratification was held fixed throughout the study. To generate an ISW of depression a second
$\tanh$
profile was used to describe the initial density field in the horizontal direction (see step in density at
$x\approx -5.2$
m) and as such mimics a ‘lock release’ experiment like those used in the laboratory (Carr et al. Reference Carr, Stastna and Davies2010; Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022). The step in the density field across the ‘lock’ was varied in vertical extent, resulting in waves of differing amplitude. The numerical domain was chosen to have dimensions that can be compared with previous experimental and numerical studies (Boegman et al. Reference Boegman, Ivey and Imberger2005; Aghsaee et al. Reference Aghsaee, Boegman and Lamb2010; Carr et al. Reference Carr, Stastna and Davies2010; Sutherland et al. Reference Sutherland, Barrett and Ivey2013; Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022, Reference Hartharn-Evans, Stastna and Carr2023; Yuan et al. Reference Yuan, Luo, Wang, Lin, Liu and Song2024), and match exactly the domain presented in Hartharn-Evans et al. (Reference Hartharn-Evans, Carr, Stastna and Davies2022), with the exception of the undulating topography.
An undulating linear topographic slope was modelled using the following profile (Lamb & Nguyen Reference Lamb and Nguyen2009):
\begin{equation} z = s \left (\hbox{itanh}(x, L_x - L_s, d) - \hbox{itanh}(x, L_x, d)\right ) +a_u\sin (f\pi x)e^{-4s_u(x-p_u)^2}, \end{equation}
where
\begin{equation} {\rm{itanh}}(x,L,d) = {1 \over 2}\left( {x - L + d\ln \left( {2\cosh \left( {{{x - L} \over d}} \right)} \right)} \right) \end{equation}
and
$d$
represents a characteristic distance for the transition from a flat bed to a constant slope,
$L_s$
is the length of the slope,
$a_u$
is the amplitude of the largest undulation,
$f$
is the frequency of the undulations,
$s_u$
is the spread of the undulations and
$p_u$
is the position of the centre of the undulations. The function smooths the transition from the flat bed to the slope and adds undulations to the slope, via the first and last terms in (2.6) respectively.
No-slip boundary conditions were applied at the flat upper and mapped lower boundaries. A mapped Chebyshev grid is employed in the vertical, implying a clustering of points near both the upper and lower boundaries that scales with the number of points in the vertical squared and that the vertical resolution improves over the slope. Free-slip boundary conditions were applied at the vertically oriented left and right ends of the computational domain, the grid spacing of which was regularly spaced.
Grid resolution was 4096 points in the
$x$
coordinate and 256 grid points in the
$z$
coordinate, giving
${\rm d}x = 1.7$
mm for simulations in a
$7$
m long tank (reduces to
${\rm d}x = 3.5$
mm in extended
$14.5$
m tank length simulations), and away from the slope,
${\rm d}z$
varies between
$0.124$
mm in the BBL and
$1.8$
mm near mid-depth. This is comparable resolution to past studies (e.g. Xu & Stastna Reference Xu and Stastna2020; Hartharn-Evans et al. Reference Hartharn-Evans, Stastna and Carr2023) and matches the resolution of past simulations shown to be in excellent agreement with corresponding laboratory experiments (Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022). Three-dimensional extensions are carried out for a subset of the cases. These are spawned from the 2-D runs prior to significant overturning. The spanwise extent of the simulations is
$0.1$
m, with 128 grid points giving a resolution of
$0.8$
mm. While comparing different types of numerical model (e.g. pseudospectral to finite volume) the 3-D simulations compare favourably with published work (e.g. Arthur & Fringer Reference Arthur and Fringer2014). The focus of the present work is not on the 3-D simulations, and more detailed analysis will be presented in future work.
2.2. Parameter variation
In the plunging, collapsing and surging cases the slope steepness and domain length were fixed at
$s=0.2$
and
$L_x=7$
m respectively. Different shoaling types were generated by varying the wave amplitude,
$a$
, (via the size of the vertical step in density over the ‘lock’) and hence wave steepness
$a/\lambda$
, where
$\lambda$
is half-wavelength as defined in Hartharn-Evans et al. (Reference Hartharn-Evans, Carr, Stastna and Davies2022). In the fission case, the slope steepness and domain length were fixed at
$s=0.033$
and
$L_x=14.5$
m respectively and wave steepness was held fixed at
$a/\lambda =0.07$
. A total of 40 model runs were carried out; 10 in each of the plunging (
$a/\lambda = 0.1, s = 0.2$
), collapsing (
$a/\lambda = 0.08, s = 0.2$
), surging (
$a/\lambda = 0.01, s = 0.2$
) and fission (
$a/\lambda = 0.07, s = 0.033$
) regimes. In all cases the stratification in the main section of the domain was held fixed at the values given above and depicted in figure 1. The main variables throughout the study were the shoaling type, undulation maximum amplitude
$a_u$
and undulation frequency
$f$
. The simulation parameters match simulations
$26\_091120$
(collapsing),
$28\_121120$
(fission),
$27\_111120$
(surging) and
$24\_071020/8L$
(plunging) presented in Hartharn-Evans et al. (Reference Hartharn-Evans, Carr, Stastna and Davies2022), where
$a_u=0$
. A smooth topographic slope is attained by simply setting
$a_u=0$
. Note that in the plunging, collapsing and surging cases the undulations on the slope typically had length scales that were four times smaller than the wavelength. In the fission case, since the slope was longer, it was possible to vary the topography parameters more widely. Consequently undulations in the fission case can have length scales that are larger than the wavelength. The distinction between these two regimes for the fission case are given in § 3.1.4. Table 1 and figure 2 provides a summary of the range of the topographic variables investigated.
Shoaling type P – plunging, C – collapsing, S – surging, F – fission and variables for the maximum undulation amplitude
$a_u$
and undulation frequency
$f$
in the undulating cases. Smooth cases for P, C, S and F had
$a_u=0=f$
.
Table 1. Long description
A table with three columns and six rows. The columns are labeled ‘Shoaling Type’, ‘a_u’, and ‘". The "Shoaling Type’ column includes entries for ‘P, C and S’ and ‘F’. The ‘a_u’ column lists various values such as ‘0.001, 0.005, 0.01, 0.015, 0.02’, ‘0.01’, ‘0.001, 0.01’, ‘0.06’, ‘0.03, 0.06’, and ‘0.01, 0.015, 0.02, 0.03’. The ‘" column lists values such as "10’, ‘1, 5, 10, 15, 20’, ‘10’, ‘0.2’, ‘0.5’, and ‘0.6’. The table compares different shoaling types and their corresponding variables for maximum undulation amplitude and frequency.", "EDH": "
Overview of topographic parameter space investigated. Shoaling type is denoted P-plunging, C-collapsing, S-surging and F-fission.
Previous work (Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022) has shown excellent agreement between the numerical model and laboratory observations in the smooth slope case (
$a_u=0$
) providing physical validation of the use of the numerical code here. Moreover, the code has been thoroughly validated in a number of other configurations, including settings that involve ISW-induced flow over irregular topography (e.g. Carr et al. Reference Carr, Stastna and Davies2010; Deepwell et al. Reference Deepwell, Stastna, Carr and Davies2017, Reference Deepwell, Stastna, Carr and Davies2019; Posada-Bedoya et al. Reference Posada-Bedoya, Olsthoorn and Boegman2025). The numerical code is available for download through its online manual: https://wiki.math.uwaterloo.ca/fluidswiki/
3. Results
3.1. Shoaling dynamics
In this section, the effect of variable topography on the different shoaling dynamics is illustrated. In subsections (3.1.1), (3.1.2), (3.1.3) and (3.1.4), the plunging, collapsing, surging and fission dynamics are, respectively, presented.
3.1.1. Plunging
In figure 3 the effect of variable topography on the plunge dynamics is shown. The colour scheme depicts vorticity,
$\omega$
, and grey contours are lines of constant density within the pycnocline. The first, second and third columns correspond to maximum undulation amplitudes of
$a_u=0$
,
$0.001$
and
$0.01$
m, respectively. When the topographic slope is smooth (figure 3
a–e), a plunge dynamics, in keeping with past studies (Sutherland et al. Reference Sutherland, Barrett and Ivey2013; Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022), is seen. Namely, the leeward face of the wave steepens, becomes gravitationally unstable and plunges over upon itself with a mushroom-like shape (figure 3
c–e). As the wave shoals, boundary-layer separation occurs in the adverse pressure gradient region, and a vortex forms within the separated flow; see growth of the blue vorticity signal at the bed, figure 3(c–e). In addition, an oppositely signed (red) vortex forms within the pycnocline due to the wave plunging dynamic. The two vortices stay distinct and the boundary-layer vortex does not affect the wave form. This is in contrast to the collapsing case (see § 3.1.2) and is a means of distinguishing between these two types of shoaling dynamics.
Sequential plots of plunging breaker types over three slopes with increasing amplitudes of undulation. The colour scheme depicts vorticity
$\omega$
and grey contours are lines of constant density within the pycnocline. The frequency of undulations is held fixed at
$f=10$
m
$^{-1}$
and in the first, second and third columns (panels (a)–(e), ( f)–( j) and (k)–(o)), the maximum undulation amplitude
$a_u=0$
,
$0.001$
and
$0.01$
m, respectively. The wave steepness is
$a/\lambda =0.1$
. The time interval between panels is
$1$
s. Supplementary movies available online.
The effect of adding undulations of varying amplitude to the slope can be seen in figures 3( f–j) (
$a_u=0.001$
m) and 3(k–o) (
$a_u=0.01$
m). In both cases, the form of the topographic bed is felt by the shoaling wave and the pycnocline is deformed accordingly. This is illustrated in the small-amplitude undulation case upon close inspection of figure 3( f–j) and can clearly be seen in the large-amplitude case figure 3(k–o). Consequently, the waves over the undulating topography do not plunge in the same clean mushroom-like fashion as that displayed in the smooth bed case (cf. figures 3
e, 3
j and 3
o) and small-scale chaotic-vortex-like features are seen that are not evident in the smooth case. In addition, overturning Kelvin–Helmholtz (K-H)-type billows can be seen on the pycnocline in the largest undulation case (figure 3
k–m) due to the interaction of the pycnocline with the large corrugated bed. The corrugations appear to set up resonant vorticity waves that cause undulations in the pycnocline and have the same wavelength as the corrugations; the K-H billows are paired with the corrugations. In figure 4, Richardson number plots confirm that the billows are highly likely a consequence of shear-induced instability. Figure 4(a) is displayed
$2$
s before figure 3(k) (i.e. prior to strong wave interaction with the bed) and (b) is displayed at the same time as figure 3(k). Values of
$Ri\in [0,0.25]$
(a condition that is necessary but not sufficient for shear instability) are given in grey scale. Finally, note that, in the largest-amplitude case, small vortices are generated at the corrugated bed behind the wave (see figure 3
k–n).
3.1.2. Collapsing
As in figure 3 but for collapsing breaker types. The wave steepness is
$a/\lambda =0.08$
. The time interval between panels is
$1$
s. Supplementary movies available online.
In figure 5 the effect of variable topography on the collapse dynamics is illustrated. When the topographic slope is smooth (figure 5 a–e), a collapse dynamics, in keeping with past studies (Boegman et al. Reference Boegman, Ivey and Imberger2005; Sutherland et al. Reference Sutherland, Barrett and Ivey2013; Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022), is seen. Namely, the leeward face steepens and begins to go gravitationally unstable (figure 5 a–c) but simultaneously the (blue) boundary-layer vortex grows in vertical extent to the point that it interacts and merges with the wave form (figure 5 c–d) and consequently pulls the wave interface back upon itself, causing the wave to collapse (figure 5 d–e).
The effect of variable topography can be seen in figures 5( f–j) and 5(k–o). Similarly to the plunging case, the form of the topographic bed is felt by the shoaling wave; the pycnocline is deformed accordingly and consequently the waves over undulating topography do not collapse in the same fashion as that displayed in the smooth bed case (cf. figures 5 e, 5 j and 5 o). Again, the larger the undulations the more disrupted the collapse dynamics is and the higher the small-scale chaotic-like motion is. Note the number of K-H billows, in the largest-amplitude case, are fewer and less well defined than in the plunging case.
3.1.3. Surging
In figure 6 the effect of variable topography on the surge dynamics is illustrated. In the surging case the wave has a very small amplitude and, consequently, the fluid dynamics is very gentle. Note that the colour bar scale and axes differ from the figures in § 3.1.1 and 3.1.2, in order to better highlight the surging dynamics. In the smooth and small topographic amplitude cases, as the wave shoals, a small bolus of fluid forms at the bed and surges upslope (see figure 6
a–e, f–j). In the largest-amplitude case (figure 6 k–o), the topography causes distortion to the density field where the pycnocline intersects the slope. Some of the incident wave energy is reflected off the slope and the amount of fluid that surges upslope is less than in the smoother cases, the topography acting like a partial barrier. This case is considered in more detail at later time steps in figure 7. The density anomaly field (
$\rho '=\rho -\rho _1$
) is displayed in (a)–(e) and the vorticity field in ( f)–( j) with grey lines of constant density superimposed. It can be seen that the reflected wave signal takes the form of a mode-2 internal wave train propagating offshore (figure 7
d,e and 7
i, j).
Sequential plots of the density anomaly
$\rho '=\rho -\rho _1$
and vorticity
$\omega$
fields for the surging breaker in the largest undulation amplitude case. Panels (a) and ( f) correspond to panel (o) in figure 6. The time interval between panels is
$4$
s. Supplementary movies available online.
3.1.4. Fission
In this section the effect of variable topography on the fission dynamics is illustrated. The slope steepness,
$s=0.033$
, is less than in the plunging, collapsing and surging cases, and the horizontal domain length
$L_x=14.5$
m is extended accordingly. Consequently, a broader range of topographic parameters can be considered and in particular undulations with length scales that exceed the wavelength are presented in addition to simulations in which the topographic parameters mimic those of previous sections. The slope steepness is more representative of the field than the cases previously discussed.
In figure 8 the same topographic parameters are used as in previous sections. In the smooth slope case (panels a–d), as the internal wave shoals, the depth ratio between the layers above and below the pycnocline shifts such that the lower layer becomes thinner than the upper. This marks a turning point in the wave’s evolution. The wave’s leeward face steepens, leading to the formation of a wave of elevation that detaches from the rear of the main depression wave and travels upslope. A strong near-bed vorticity signal develops within the core of this elevation wave. This dynamics is consistent with previous findings (e.g. Xu & Stastna Reference Xu and Stastna2020; Ghassemi et al. Reference Ghassemi, Zahedi and Boegman2022; Hartharn-Evans et al. Reference Hartharn-Evans, Carr, Stastna and Davies2022, Reference Hartharn-Evans, Stastna and Carr2023). As in the plunging and collapsing cases, when the topography is of sufficient amplitude the wave itself appears to be deformed by interaction with the undulating bed enhancing vorticity across the pycnocline prior to any BBL-derived vortices shedding upwards and interacting with the pycnocline (panel i). In the largest-amplitude case, this drives overturning K-H billow-type structures (panels i–k). In addition there is interaction of BBL-derived vortices with the pycnocline, which leads to enhanced small-scale chaotic-like features in the wave tail.
As in figure 3 but for fission breaker types. The wave steepness is
$a/\lambda =0.07$
. The time interval between panels is
$8$
s. Supplementary movies available online.
In figure 9 the frequency of the undulations is reduced to
$f=0.5$
m
$^{-1}$
and the maximum amplitude is
$a_u=0.03$
m. In this case the length scale of the topography exceeds that of the wavelength. In the variable topography case (panels f–g), initially, a wave of elevation begins to form similarly to the smooth slope case (panels a–e). However, as the wave moves into regions of locally deeper water – where the topographic length scale exceeds that of the wave–the lower layer no longer remains thinner than the upper layer. This stops the classical fission process observed for smooth topography, and instead the fluid adjusts by generating an oscillatory tail behind the main depression wave (panels g–j). Additionally, finer-scale motion develops along the pycnocline in the wave’s wake, and wave energy is scattered at the topographic bump via small-scale chaotic-like features. Notably, the coherent BBL elevation wave dynamics seen in the smooth case (panel e) is absent in the undulating case (panel j). The shoaling dynamics is locally sensitive to the topographic form (i.e. non-monotonic slopes), as noted for free-slip simulations on the field scale by Lamb & Warn-Varnas (Reference Lamb and Warn-Varnas2015).
Sequential plots of fission breaker types over a smooth (panels a–e) and undulating (panels f–j) slope. The frequency of the undulations is
$f=0.5$
m
$^{-1}$
and the maximum undulation amplitude is
$a_u=0.03\,\text{m}$
. The wave steepness is
$a/\lambda =0.07$
. The time interval between panels is
$3$
s. Supplementary movies available online.
Figure 10 presents the later evolution of the undulating slope case. The first panel corresponds to the final panel of figure 9, with the
$x$
-axis shifted to the right. As the wave continues to propagate and the local water depth changes, once more–passing the turning point, a second fission event occurs (see panels e and f). The fact that the passage through the turning point is arrested over a significant portion of the domain prior to this leads to the generation of a significant wave train. In the terminal phase of shoaling, when the pycnocline intersects the slope, the wave of elevation transitions into boluses of fluid that travel upslope. This behaviour aligns with the terminal shoaling dynamics observed in the smooth case, as described by Hartharn-Evans et al. (Reference Hartharn-Evans, Stastna and Carr2023) and Castro-Folker & Stastna (Reference Castro-Folker and Stastna2024).
4. Three-dimensional effects
4.1. Methodology
In order to assess (i) the ability of 2-D simulations in capturing the onset of breaking and (ii) the effect of undulating topography on three-dimensionalisation, a limited number of 3-D simulations were performed. Three-dimensional simulations were carried out in the plunging regime with undulating topography and the collapsing regime with and without undulations. These simulations were spawned from the 2-D simulations
$1$
s earlier than the earliest panels shown in figures 3 and 5, respectively. The purpose of the 3-D simulations was to ensure that the qualitative aspects of the 2-D simulations persist in the 3-D situation. This proved to be the case for all three cases carried out. In particular, the upslope edge of the partially broken down wave remained remarkably consistent between 2-D and 3-D simulations even after significant 3-D effects set in downslope.
4.2. Results
Figure 11 summarises the results for the collapsing case. Panels (a) and (b) show the shaded distribution of the mean of
$v^2$
, where
$v$
is the spanwise velocity, saturated at
$50$
$\%$
of the maximum. The primary effect of the topography is to yield a greater along-slope extent of the spanwise component of the velocity. Panel (c) depicts the 3-D wave form in a subdomain over the slope. The isopycnal at the centre of the pycnocline is shown in red, and the undulating slope is shown in yellow. Upslope of the breaking region the topography is masked by the red isosurface since the wave-induced currents draw down the pycnocline so it undulates in close proximity of the topography. Small-scale structure is associated with the downstream face of the wave, and extends downslope from this region. There is a clear separation of scales between the scale of the visible 3-D variations and the larger scale associated with the wave. The large-scale structure of the wave tracks the 2-D simulations surprisingly well for such a mature, three-dimensionalised state.
Summary figure of the 3-D collapsing simulations. The mean of
$v^2$
(where
$v$
is the spanwise velocity) scaled by the maximum value of the smooth case for the (a) smooth slope case and (b) undulating case (
$a_u=0.01$
m). (c) An isosurface of density (red) for the undulating topography (yellow) case. The value of density chosen is that for the centre of the pycnocline. (d) Maximum values of
$|v|$
for the cases with undulating topography (blue) and smooth slope (red), the time of panels (a–c) is indicated by a vertical, dashed black line.
Figure 11(d) shows the time series of the maximum of
$|v|$
for the case with topography (blue) and the smooth slope case (red). It can be seen that the effect of the topography is to induce faster three-dimensionalisation. The maximum of
$|v|$
for the smooth slope case is still increasing at the final time shown, while the case with topography peaks at the time shown in panels (a–c), which is indicated by a vertical, dashed black line in (d).
The plunging case (not shown) yields a much larger breaking region. The 3-D run remains consistent in large-scale structure with the 2-D case up to the final panel shown in figure 3. The primary way in which small-scale vortices are modified from the 2-D simulation is by spilling downslope beyond the downstream breaking face of the wave, consistent with Xu et al. (Reference Xu, Stastna and Deepwell2019) and the terminal or bolus stage of Castro-Folker & Stastna (Reference Castro-Folker and Stastna2024).
Figure 12 returns to the collapsing case and shows the along-slope (a,b) and spanwise (c,d) components of the bottom shear. The smooth case is shown in the left column (a,c) and the case with undulating topography is shown in the right column (b,d). The subset of the domain that is shown corresponds to that shown in figure 11(a–c) at the same time. The colour bar is saturated at
$\pm 25$
$\%$
of the overall maximum value of the shear. The along-slope component of shear has a clear footprint of the main wave, and in the case with undulating topography, modulations due to the individual topographic hills and valleys (b). It can be seen that the 3-D structure is more developed in panels (b,d), and indeed the spanwise component of the shear for the smooth slope case (panel c) is quite weak. For the case with undulating topography the along-slope component of shear (b) has some medium-scale structure (3–4 peaks) in the spanwise direction near
$x=0.82$
m. The spanwise component of shear (d) is ‘splotchy’ in nature with both positive and negative regions visible. The fact that the magnitude of the spanwise shear is similar to that of the along-slope shear suggests that any sediment resuspension parametrisations should account for both components of bottom shear.
The components of the bottom shear at the time denoted by the vertical, dashed black line in figure 11 (d). The along-slope (a,b) and spanwise (c,d) components of the bottom shear for the case with undulating topography (a,c) and without (b,d) are shown, saturated at
$25$
$\%$
of the maximum absolute value of the shear.
For the times discussed there was no evidence of the coherent spanwise structures reported in the measurements off the Scripps Pier by Lucas & Pinkel (2022). Extending the spanwise extent of the 3-D simulations is a clear direction for future work. Similarly, a detailed analysis of the fully three-dimensionalised state is left for a future publication.
5. Potential regions of mixing identified via the USP method
Typically in process-based studies, energetic differences between simulations can be analysed using measures such as the bulk viscous dissipation rate, and the diapycnal mixing rate. The most common way of discussing the latter requires comparison of density profiles with a ‘sorted’ density profile to split the potential energy into available potential energy and background potential energy (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995). The sorted profile is unambiguously defined when there is no bathymetry, but complications arise for complex bathymetry, making this approach unsuitable for this setting (Hartharn-Evans et al. Reference Hartharn-Evans, Stastna and Carr2024). Moreover, local information about where mixing occurs, and what dynamical features it is associated with, is not immediately clear in the sorting formalism. The USP method of Hartharn-Evans et al. (Reference Hartharn-Evans, Stastna and Carr2024) allows for the identification of regions where various scalars (e.g. density, kinetic energy, enstrophy) attain user-defined properties. These can then be tracked back to physical space, to identify dynamic features that indicate a high likelihood of mixing. In the following, two-variable probability density histograms of the magnitude of the density gradient
$|\boldsymbol{\nabla }\rho |$
and a kinematic tracer were used to identify fluid regions of interest (ROIs) which have particular fluid properties, rather than being geometrically defined. The kinematic tracer chosen was enstrophy,
${{\varOmega}} = ({1}/{2})|\nabla \times {\bf{u}}{|^2}$
, a quantity that describes rotational energy, and which for 2-D flows is related to dissipation. As such, it is associated with (but is not a direct measure of) diapycnal mixing. The exact choices of parameters for ROI plots (and possible alternatives), and the underlying changes in paired histogram plots are given in further detail in Appendix A, but in summary are chosen to represent paired fluid properties that are not observed prior to the wave breaking dynamics.
Plots indicating regions of potential mixing for three indicative time steps (
$t = 52, 59, 64$
s respectively) during the flow evolution of the plunging breaker over a smooth (
$a_u = 0$
), small-amplitude (
$a_u = 0.001$
m) and large-amplitude (
$a_u=0.01$
m) undulating topography. The physical location of potential mixing regions (identified from ROI double histogram plots) are displayed in blue on top of the enstrophy field (orange colour bar). The ROI is defined as
$|\boldsymbol{\nabla }\rho | \geqslant 3000$
kg m
$^{-4}$
and
$\varOmega \geqslant 58$
s
$^{-2}$
.
Regions of potential mixing as in figure 13, but for the collapsing case. The ROI is defined as
$|\boldsymbol{\nabla }\rho | \geqslant 3000$
kg m
$^{-4}$
and
$\varOmega \geqslant 26$
s
$^{-2}$
.
Plots showing the spatial location of fluid for a given ROI are given in figures 13, 14 and 15 in the plunging, collapsing and fission regimes, respectively. The location of the fluid in physical space identified from the ROI is displayed in blue on top of an orange colour map marking the enstrophy for context. In the first row of all plots, the wave is displayed prior to interaction with the bed and no fluid is identified as being within the ROI at this stage. As such the first row of the plots in physical space do not contain any ‘blue’ fluid and can be thought of as the base case prior to any potential mixing events. Across all cases of surging, very few regions where mixing may occur were identified, and these plots are not shown for brevity.
Regions of potential mixing as in figure 13, but for the fission case and for the smooth and large
$a_u$
cases only. The indicative time steps in this case are
$t = 80, 104, 124$
s respectively. The ROI is defined as
$|\boldsymbol{\nabla }\rho | \geqslant 3000$
kg m
$^{-4}$
and
$\varOmega \geqslant 18$
s
$^{-2}$
.
The USP analysis reveals extended regions where mixing may occur (henceforth referred to simply as mixing) over undulating slopes compared with the smooth slope during the early stages of plunging breaking (figure 13 d–f), especially for the largest undulation. This enhancement is due to distortion of the pycnocline (due to K-H shear instability, as seen in § 3.1.1) above each topographic bump (figure 13 f). While the wave is breaking (figure 13 d–f), the differences in the mixing dynamics are clear from the differences in regions of mixing. In the smooth and small-amplitude undulation cases, as the pycnocline plunges large coherent regions of mixing are confined to the pycnocline and follow the mushroom-like shape associated with plunging (figure 13 d–e). This is represented broadly by enhanced enstrophy at the pycnocline. In contrast, the large-amplitude undulation induces multiple small billows on the pycnocline very close to the boundary (figure 13 f), and as such mixing is concentrated in two regions: near the boundary layer, as well as higher in the water column, with small-scale patches of mixing expected throughout the region occupied by the breaking wave. Nonetheless, late stage mixing characteristics, after the wave has passed (figure 13 g–i), are remarkably similar across all roughness cases.
The regions of potential mixing for collapsing are identified in figure 14 and show similar trends to the plunging case. Namely, during the collapse stage (figure 14 d–f) the smooth and small-amplitude cases are very similar, with mixing regions primarily due to the collapse dynamics, while in the undulating amplitude case, the mixing region exhibits differences due to the formation of K-H billows within the pycnocline. At later times, after the wave has passed (figure 14 g–i) there is little difference between the three cases.
Figures 15 presents the USP plots for the fission case. The fission case shows less/smaller regions of potential mixing (blue fluid) than the plunging and collapsing cases, which is to be expected as the fission case is less dynamic (cf. results in §§ 3.1.1, 3.1.2 and 3.1.4). During the fission process (panels (d) and ( f)), the highest-amplitude undulating case exhibits enhanced enstrophy at the bed (due to vortices formed within the corrugations) and across the pycnocline (due to K-H instability) compared with the smooth case.
6. Conclusion
It has been shown that variable topography can fundamentally change the shoaling dynamics in all four breaking cases; plunging, collapsing, surging and fission. Notably in the plunging, collapsing and fission cases, if the corrugations at the bed are large enough, they can deform the wave shape and induce K-H billows on the pycnocline. In the surging case, undulating topography can reduce upslope surging and cause a mode-2 wave signal to be reflected offshore. The USP analysis showed that potential local mixing processes are affected by variable topography, but the differences are transient and restricted to the breaking phase of shoaling, with long-term effects on possible mixing regions surprisingly unaffected by variable topography.
The presence of an undulating topography modifies the local water layer depths which are a fundamental factor in the fission process, and as such the fission dynamics can be significantly modified in response if the length scale of the topography exceeds that of the wavelength. In this case, local changes in topography can arrest the passage through the turning point. In this situation, a new wave train develops over the deeper region and, as its crests shoal, they may induce small-scale chaotic-like features and potential mixing at a later time and in a manner in which inter-crest interaction may become important (see figure 10 i–j). In the field, where fission dominates shoaling processes, topographic features are expected to locally modify the shoaling process and, if large enough, could effectively block the transport of energy and scalar properties further upslope, with implications for the upslope environment.
Recent laboratory work by Wang et al. (Reference Wang, Yuan, Lin and Song2025) on a reef–lagoon system has concentrated on steep slopes (
$s=0.5$
and
$s=0.2$
) with roughness elements added to mimic coral. The flow induced far more near bottom turbulence than the simulations presented here. A detailed inter-comparison with the present work is a clear avenue for future work, although the detailed dynamics in the lagoon is not shown in Wang et al. (Reference Wang, Yuan, Lin and Song2025).
Although this manuscript focusses on 2-D simulations, 3-D simulations (§ 4) and Castro-Folker & Stastna (Reference Castro-Folker and Stastna2024) suggest that the overall dynamics remains qualitatively consistent. In particular, Castro-Folker & Stastna (Reference Castro-Folker and Stastna2024) found that boundary-layer dynamics during fissioning is dominated by 2-D separation, with significant three-dimensionalisation occurring only when ejected vortices interact with the pycnocline downslope of the wave and in the terminal bolus stage of shoaling. Papoutsellis et al. (Reference Papoutsellis, Lucas, Shukla, Okun, Pinkel and Schmidt2026) carried out 3-D field-scale simulations of obliquely shoaling ISWs using the MITgcm model. Their simulations showed upslope propagating boluses of cold water during the fission process exhibiting striking qualitative agreement with the experimental observations of Ghassemi et al. (Reference Ghassemi, Zahedi and Boegman2022), and the 2-D direct numerical simulation results in Hartharn-Evans et al. (Reference Hartharn-Evans, Stastna and Carr2023). For the present simulations the along-slope distribution of the spanwise average of
$v^2$
and the 2-D enstrophy USP plots were compared. It was found that the distribution in the two plots coincides to a significant degree even when the rear of the wave has fully three-dimensionalised. The USP plots show unphysical local patches of enstrophy that are smeared out in the 3-D simulations. However, as a means of identifying where mixing will most likely take place in 3-D simulations the USP plot provides a very useful tool.
Field conditions are often far more complex, with significant background turbulence. Zulberti et al. (Reference Zulberti, Jones and Ivey2020) reported modulation of background turbulence by mode-1 ISWs on the Northwest Australian shelf, finding that the wave increased turbulent boundary-layer height near the wave crest. High levels of background turbulence in such situations are difficult to replicate in laboratory settings. Posada-Bedoya et al. (Reference Posada-Bedoya, Olsthoorn and Boegman2024) proposed mimicking the background turbulence in the boundary layer by having an ISW repeatedly pass through its own wake and boundary-layer instability that slowly develops within it. The repeated passage shifted the locus of instability closer to the wave crest. It is unclear how to carry these results over to shoaling, where repeated passage is impossible.
In subsequent work, Posada-Bedoya et al. (Reference Posada-Bedoya, Olsthoorn and Boegman2025) modified the instability by introducing small-scale undulations (albeit two-dimensional) mimicking laboratory conditions. The undulations seeded noise and the ISW acted as an amplifier for certain portions of the noise spectrum. Introducing sandpaper-scale roughness in addition to the topographic features studied here is an interesting avenue for future work, although it should likely be pursued with a numerical model tailored to fully 3-D roughness.
In measurements off the Scripps Pier, Lucas & Pinkel (2022) and Papoutsellis et al. (Reference Papoutsellis, Lucas, Shukla, Okun, Pinkel and Schmidt2026) reported coherent pulses during shoaling (after passage through the turning point in the fission case). Determining whether coherent structures in the spanwise direction are possible in the collapsing or plunging cases would require wider domains than those used herein, and is a clear avenue for future work.
Scaling the results presented here to the field requires caution. For plunging and collapsing cases, turbulence production is expected to be the most important consideration and the insensitivity of the final state to bathymetric details is encouraging. The fissioning case, which is likely more relevant to most field situations, is more subtle. There is existing evidence that large-scale, non-monotonic features (Lamb & Warn-Varnas Reference Lamb and Warn-Varnas2015) can fundamentally change wave evolution during shoaling before boundary-layer effects begin to be significant, and this study adds evidence for the portion of shoaling in which the boundary layer is active.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.12129.
Acknowledgements
The authors acknowledge the role played by the regular series of Scotland-Norway Waves symposia in helping to facilitate the collaborative work reported in this paper and providing a platform for disseminating the results. The authors thank Prof Pete Diamessis and two anonymous reviewers for their critical reviews which have led to improvements in the paper.
Funding
The School of Mathematics, Statistics and Physics at Newcastle University are thanked for supporting Jack Smith through a summer studentship. Samuel Hartharn-Evans was supported by the Natural Environment Research Council (NERC) funded ONE Planet Doctoral Training Partnership (grant number NE/S007512/1) and a Turing Global Fellowship. This work made use of the Rocket High Performance Computing service at Newcastle University. Partial funding for travel was provided by the Natural Science and Engineering Research Council of Canada RGPIN-311844-37157.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available on Newcastle Universities data repository: https://doi.org/10.25405/data.ncl.30847469.
Appendix A. The USP methodology
The USP method introduced in Hartharn-Evans et al. (Reference Hartharn-Evans, Stastna and Carr2024) provides a qualitative method to identify regions where mixing may take place, as an alternative to quantitative methods that may not be valid for all scenarios e.g. when complex topography is involved. The following appendix provides a worked example based on the plunging simulation in § 3.1.1 where the undulation amplitude
$a_u=0.01$
m.
First, we choose a pair of variables to consider as combined proxies for mixing. The first is typically a tracer-based measure (e.g.
$\rho$
,
$|\boldsymbol{\nabla }\rho |$
,
$\partial \rho /\partial z$
, temperature or salinity), whilst the second variable is a kinetic measure (e.g. kinetic energy density, enstrophy, dissipation). The user could also consider implementing other variables that can be linked to the energetics and physics of mixing, especially for fully resolved 3-D simulations, although these introduce additional complexities. There are arguments for a range of different measures that the user can weigh up on a case-by-case basis:
-
(i) Density alone is a very simple measure, and changes in density distributions are the thing we are actually interested in above almost anything else in a mixing study.
-
(ii) Vertical density gradient brings in some elements of the total gradient of density, whilst highlighting regions of overturning that are inherently linked to mixing. In conjunction with a relevant kinematic variable (e.g. velocity shear), it can give us insights directly into the Richardson number-type dynamics.
-
(iii) Density gradients consider each direction and the total action of stirring.
-
(iv) The Laplacian of the density would give the part of the heat equation that diffusivity acts directly upon.
In this study,
$|\boldsymbol{\nabla }\rho |$
, is chosen as the tracer-based measure. It is a representation of the tracer gradients that diffusion acts on. Enstrophy,
$\varOmega$
, is chosen as the kinematic measure, it represents the kinetic action of strain structures. Taken together, these variables might tell us where the action of the kinetics is straining density gradients and potentially enhancing regions of mixing.
Snapshots of USP variables
$|\boldsymbol{\nabla }\rho |$
and
$\varOmega$
at two time steps relevant to potential mixing regions for the plunging simulation with undulation amplitude of 0.01m, prior to interaction with the slope (left-hand side,
$t = 40$
s), and once breaking is underway (right-hand side,
$t = 60$
s).
Inspection of the raw data fields helps us identify relevant spatial limits and variable limits. Many variables have some very long tails on their distributions, and so cutoff at a reasonable level is required. We show two time points in figure 16, the first is just prior to the wave interacting with the slope (as close to fully developed flow prior to the features we are interested in), and a later one chosen roughly as the wave interacts with the slope.
-
(i) In the early plot,
$|\boldsymbol{\nabla }\rho |$
is a signifier of the pycnocline, and we can see an area of small-scale vortices in the wave generation region. -
(ii) At this time, enstrophy is high in a number of locations, most prominently it outlines the wave itself, notably not throughout the entire waveform of the pycnocline. It also outlines the tail of the wave, but not ahead of it (unlike
$|\boldsymbol{\nabla }\rho |$
). Finally, there are regions at the upper and lower no-slip boundaries above and below the wave itself that do not align with similarly high
$|\boldsymbol{\nabla }\rho |$
regions. -
(iii) At later times,
$|\boldsymbol{\nabla }\rho |$
is enhanced in what we identify as the wave breaking region, reaching higher values than we saw at earlier time steps. -
(iv) Similarly, enstrophy is broadly enhanced in the wave breaking region.
From this point, we can begin to consider the paired histogram elements of this method. What might we expect based on the above plots?
-
(i) The majority of the volume of fluid is represented by both
$|\boldsymbol{\nabla }\rho |\approx 0$
and
$\varOmega \approx 0$
. -
(ii) A region of elevated enstrophy with
$|\boldsymbol{\nabla }\rho |\approx 0$
is associated with the boundary layers. -
(iii) A region of elevated
$|\boldsymbol{\nabla }\rho |$
with
$\varOmega = 0$
is associated with the pycnocline at rest, prior to the wave passing. -
(iv) Finally, a region where both enstrophy and
$|\boldsymbol{\nabla }\rho |$
are elevated. This is the bit of the regime that is most interesting from a mixing perspective.
Paired histograms showing the probability density function of the combination of two variables at two time steps prior to mixing (a) and once breaking is underway (b). White boxes outline the relevant ROI for this flow, identified from the extremes in (a).
We are interested in what mixing occurs beyond the background level (as shown in figure 17). With this choice of variables, we know that ‘mixing’ would be associated with elevated enstrophy, and elevated
$|\boldsymbol{\nabla }\rho |$
. Currently, the USP method requires a rectangular USP domain for the ROI, based on the above plot, we might reasonably choose
$|\boldsymbol{\nabla }\rho |\geqslant 3000$
kg m
$^{-4}$
and
$\varOmega \geqslant 50$
s
$^{-2}$
.
Once the ROI is confirmed, we can move back to the physical plots, and identify where in the flow the fluid represented by the ROI in USP space is in physical space (indicated by blue in figure 13 for example).
ρ′=ρ−ρ1
ρ1
ρ
au
f
au=0=f
ω
f=10
−1
au=0
0.001
0.01
a/λ=0.1
1
Ri
2
a/λ=0.08
1
a/λ=0.014
3
ρ′=ρ−ρ1
ω
4
a/λ=0.07
8
f=0.5
−1
au=0.03m
a/λ=0.07
3
5
v2
v
au=0.01
|v|
25
%
t=52,59,64
au=0
au=0.001
au=0.01
|∇ρ|⩾3000
−4
Ω⩾58
−2
|∇ρ|⩾3000
−4
Ω⩾26
−2
au
t=80,104,124
|∇ρ|⩾3000
−4
Ω⩾18
−2
|∇ρ|
Ω
t=40
t=60
Open access
