Surface wavepackets subject to an abrupt depth change. Part 1. Second-order theory

Abstract

This paper develops second-order theory for narrow-banded surface gravity wavepackets experiencing a sudden depth transition based on a Stokes and multiple-scales expansion. As a wavepacket travels over a sudden depth transition, additional wavepackets are generated that propagate freely obeying the linear dispersion relation and arise at both first and second order in wave steepness in a Stokes expansion. In the region near the top of the depth transition, the resulting transient processes play a crucial role. At second order in wave steepness, free and bound waves coexist with different phases. Their different speeds of travel result in a local peak a certain distance after the depth transition. This distance depends on the water depth $h_s$ relative to the carrier wavelength on the shallower side $\lambda _{0s}$. We validate our theory through comparison with fully nonlinear numerical simulations. Experimental validation is provided in a companion paper (Li et al, J. Fluid Mech., 2021, 915, A72). We conjecture that the combination of the local transient peak at second order and the magnitude of the linear free waves provides the explanation for the rogue waves observed after a sudden depth transition reported in a significant number of papers and reviewed in Trulsen etal (J. Fluid Mech., vol. 882, 2020, R2).

1. Introduction

There is considerable engineering and scientific interest in the probability of large waves occurring in the ocean. Large waves that occur more frequently than predicted by standard linear theories are sometimes termed ‘rogue’ or ‘freak’ waves. Various physical mechanism are known to generate abnormal wave statistics as reviewed by Dysthe, Krogstad & Müller (2008), Onorato et al. (2013) and Adcock & Taylor (2014). A convenient, and commonly used, proxy for the number of rogue waves is the kurtosis (or excess kurtosis) of the free surface (Mori & Janssen 2006).

In the last decade, a number of studies have suggested that a transition of water depth could play an important role in an enhanced occurrence probability of extreme waves (Sergeeva, Pelinovsky & Talipova 2011; Onorato & Suret 2016; Trulsen 2018; Majda, Moore & Qi 2019). This phenomenon has been demonstrated both numerically (Sergeeva et al. 2011; Gramstad et al. 2013; Viotti & Dias 2014; Ducrozet & Gouin 2017; Zhang et al. 2019) and experimentally (Trulsen, Zeng & Gramstad 2012; Bolles, Speer & Moore 2019; Zhang et al. 2019; Trulsen et al. 2020). To date, a number of accidents have been reported that were seemingly caused by rogue waves in finite and shallow water depth (Chien, Kao & Chuang 2002; Nikolkina & Didenkulova 2011). This also suggests the role of a varying bathymetry in causing extreme wave events in the real world.

The mechanism causing the enhanced kurtosis at the top of slopes remains an open question, although a number of authors have pointed to the role of second-order components in wave steepness (Gramstad et al. 2013; Zhang et al. 2019; Zheng et al. 2020). Waves will interact with slopes in various ways (see e.g. Dingemans 1997; Madsen, Sørensen & Schäffer 1997; Booij, Ris & Holthuijsen 1999; Madsen & Schäffer 1999; Holthuijsen 2010). Of particular relevance for the present study are the investigations of the interplay of bound and free waves (Foda & Mei 1981; Mei & Benmoussa 1984; Battjes et al. 2004).

A useful limiting case for wave–bathymetry interaction is that of waves passing over a step, where the depth changes from one limiting (non-zero) value to a second limiting (non-zero) value, as defined in Newman (1965) and shown in figure 1. Most relevant studies have only considered linear waves and proposed various methods to deal with the presence of a step in a potential flow, as challenges exist due to the discontinuity caused by the step. Some example methods are the Green's function method proposed in Rhee (1997), wavemaker theory (Newman 1965; Havelock 1929), the long-wave approximation (Mei, Stiassnie & Yue 1989), the Galerkin-eigenfunction method (e.g. Fletcher 1984; Massel 1983, 1993; Belibassakis & Athanassoulis 2002, 2011) and direct numerical computations (Mei & Black 1969; Kirby & Dalrymple 1983). These investigations of the leading-order physics show that when the wave ‘feels’ a step in the seabed, then the wave will be partially reflected and partially transmitted. Moreover, the transmitted wave amplitude can be as large as double that of the incident wave and as small as zero in the limit in which a step becomes a wall throughout the water column (Kreisel 1949).

Figure 1. Diagram of the bathymetry and coordinate system adopted. In the diagram, we have used a narrow-banded wavepacket with surface elevation $\zeta (x) = A\exp [-{(x-x_0)^2}/{(2\sigma ^2)}] \cos (k_0x)$, where $4\sigma$ denotes the characteristic length of a Gaussian packet, $x_0=2\sigma$, $A$ is the amplitude and $k_0$ denotes the carrier wavenumber ($\lambda _0 = 2{\rm \pi} /k_0$ the carrier wavelength); $h_d$ and $h_s$ denote the water depth on the deeper and shallower sides, respectively.

For steeper waves passing over steps, second-order effects in wave steepness become significant. Massel (1983) derived second-order results for monochromatic waves. Specifically, he found that second-order superharmonic free waves are released as a result of weakly nonlinear waves interacting with a step, and the interplay of the superharmonic free and the superharmonic bound wave may result in beating near the top of the depth transition. The beating length is $2{\rm \pi} /(k_{20}-2k_0)$, in which $k_0$ is the wavenumber of the linear monochromatic wave and $k_{20}$ that of the free second-order superharmonic component, and this beating leads to a maximum of the superharmonic wave crest up to twice as large as that of the superharmonic bound wave. This beating phenomenon has been confirmed experimentally (Monsalve Gutiérrez 2017).

The present paper and its companion paper (Li et al. 2021) extend the work of Massel (1983) with the objective of explaining the mechanism behind increases in excess kurtosis observed at the top of slopes. In order to do so, this paper develops analytical solutions for narrow-banded wavepackets experiencing a sudden depth transition in the form of a step using a Stokes expansion up to second order in wave steepness. These solutions, which extend the results of Massel (1983) for monochromatic waves to wavepackets, capture the release of both sub- and superharmonic second-order free waves at the step. We validate these solutions by comparing to a fully nonlinear potential-flow model in the present paper and to experiments in the companion paper Li et al. (2021).

2. Theoretical model

2.1. Problem definition

We consider a unidirectional surface gravity wavepacket propagating in a region with an abrupt change of water depth in the framework of two-dimensional potential-flow theory, neglecting the effects of viscosity and surface tension. The bathymetry is illustrated in figure 1. The water depth $h(x)$ changes abruptly from a constant $h_d$ to $h_s$ at $x = 0$, with $x$ the horizontal coordinate. We assume $h_d \geqslant h_s$, and the water depths can be deep ($k h\gg 1$, with $k$ the wavenumber), intermediate ($k h={O}(1)$) or shallow ($k h\ll 1$) compared to the characteristic wavelength. The undisturbed water surface is located at $z = 0$. The system can be described as a boundary value problem governed by the Laplace equation:

(2.1)\begin{equation} \nabla^2 \varPhi =0\quad \textrm{for}\ -h(x)\leqslant z\leqslant \zeta(x,t), \end{equation}

where $\varPhi (x,z,t)$ is the velocity potential and $\zeta (x,t)$ is the free-surface elevation. Equation (2.1) should be solved subject to nonlinear kinematic and dynamic boundary conditions at the free surface:

(2.2a,b)\begin{equation} \dfrac{\text{D}}{\text{D}t} (z-\zeta) = 0\quad \textrm{and}\quad g\zeta+\dfrac{\partial \varPhi}{\partial t} + \dfrac{1}{2} \left(\boldsymbol{\nabla}\varPhi \right)^2=0 \quad \textrm{for}\ z= \zeta(x,t), \end{equation}

where $g$ is the gravitational acceleration; a bottom boundary condition

(2.3)\begin{equation} \dfrac{\partial \varPhi}{\partial z} =0\quad \textrm{for}\ z={-}h(x); \end{equation}

continuity of the potential and its horizontal derivative in the fluid exactly above the step

(2.4a,b)\begin{equation} \left[\varPhi\right]_{x\to 0^-} =\left[\varPhi\right]_{x\to 0^+}\quad \textrm{and}\quad \left[\dfrac{\partial \varPhi}{\partial x}\right]_{x\to 0^-} = \left[\dfrac{\partial \varPhi}{\partial x}\right]_{x\to 0^+}\ \textrm{for}\ -h_s\leqslant z \leqslant \zeta(x,t); \end{equation}

and a no-flow boundary condition on the step wall

(2.5)\begin{equation} \left[\dfrac{\partial \varPhi}{\partial x}\right]_{x\to 0^-} = 0\quad \textrm{for}\ -h_d \leqslant z<{-}h_s. \end{equation}

2.2. Stokes and multiple-scales expansions

In order to solve the boundary value problem (2.1)–(2.5), the unknown $\varPhi$ and $\zeta$ are expressed as series solutions in the wave steepness $\epsilon = k_0A$ (a so-called Stokes expansion), with $k_0$ and $A$ denoting the characteristic wavenumber and wave amplitude, respectively:

(2.6a,b)\begin{equation} \varPhi = \epsilon \varPhi^{(1)} + \epsilon^2 \varPhi^{(2)}+ {O}(\epsilon^3)\quad \textrm{and}\quad \zeta = \epsilon \zeta^{(1)} + \epsilon^2 \zeta^{(2)}+ {O}(\epsilon^3), \end{equation}

where we consider up to the first two orders. Substituting (2.6) into the the boundary value problem (2.1)–(2.5) leads to a collection of terms at the first two orders in $\epsilon$, which can be solved successively, as presented in §§ 2.5 and 2.6, respectively.

We consider a narrow-bandwidth or quasi-monochromatic wavepacket that, at least in the absence of the step, can be considered as a carrier wave whose amplitude varies slowly in both space and time (e.g. Mei et al. 1989). Both slow and fast scales are introduced in a multiple-scales expansion. Let $\psi _0 = k_0x_0-\omega _0t_0+\mu _0$ be the phase of the carrier wave, where $\omega _0$ is the angular wave frequency, $\mu _0$ is an arbitrary phase shift and $x_0$ and $t_0$ are the fast scales. We allow for slow variation of the carrier wave amplitude packet in the form of $A(X, T)$, in which $X=\delta x_0$ and $T=\delta t_0$ are the slow scales and $\delta$ is the scale separation parameter of the problem and a measure of the bandwidth of the wavepacket. In previous work, notably in Mei et al. (1989), Yuen & Lake (1975) and Dysthe (1979), the two small parameters are commonly set to be of the same order (i.e. ${O}(\delta )={O}(\epsilon )$), resulting in the derivation of a third-order packet equation of the nonlinear Schrödinger type. Herein, we do not make this assumption and focus only on the first two orders in the steepness $\epsilon$. Consequently, all the components will evolve according to the linear dispersion relationship or, for second-order bound waves, that of their linear parent waves. Derivative operators can be written in terms of a combination of fast and slow derivatives:

(2.7a,b)\begin{equation} \partial_x = \partial_{x_0} +\delta \partial_X +{O}(\delta^2)\quad \textrm{and}\quad \partial_t = \partial_{t_0} +\delta \partial_T +{O}(\delta^2). \end{equation}

Our assumption of a narrow-banded or quasi-monochromatic wavepacket that evolves slowly in time applies to the incoming and, consequently, to the transmitted and reflected wavepackets. Although the incoming, transmitted and reflected wavepackets are slowly varying in space away from the step, they are discontinuous at this location and need to be matched according to (2.4) and (2.5) to ensure continuity, resulting in the generation of evanescent waves. We examine this further below.

2.3. Description of the incoming wavepacket

Following Mei et al. (1989) and Massel (1983), we express the incoming wavepacket to leading order as

(2.8)\begin{align} \varPhi_{I} = \underbrace{\epsilon(\varPhi_I^{(11,0)} + \delta \varPhi_I^{(11,1)}+{O}(\epsilon\delta^2))}_{\epsilon\varPhi^{(1)}} +\underbrace{\epsilon^2(\varPhi_I^{(22,0)}+\delta\varPhi_I^{(20,1)} +{O}(\epsilon^2\delta^2))}_{\epsilon^2\varPhi^{(2)}} +\,{O}(\epsilon^3), \end{align}

which is valid for $x\leqslant 0$, i.e. over the flat seabed to the left of the step. The superscript $(mn,j)$ denotes the term of ${O}(\epsilon ^m\delta ^j)$ that is proportional to the harmonic $\exp (\mathrm {i} n\psi _0)$, with $n=0$ corresponding to the bound subharmonic or ‘mean flow’ and $n=2$ to the bound superharmonic (only the real part of $\exp (\mathrm {i} n\psi _0)$ is understood). An analogous equation to (2.8) describes the free-surface elevation of the incoming wavepacket $\zeta _{I}$, and we proceed to express all the solutions in terms of the packet of its lowest-order term. Specifically, we assume

(2.9)\begin{equation} \epsilon \zeta_I^{(11,0)}=A_I(X-c_{g0}T)\cos\psi_0, \end{equation}

where the amplitude packet $A_I$ is real, $c_{g0}$ is the group velocity and the dependence of $A_I$ on $X-c_{g0}T$ is based on the solvability condition (13.2.29) in Mei et al. (1989). Hence, the potentials of the incoming wavepacket at different orders are expressed as (Massel 1983; Mei et al. 1989; Calvert et al. 2019)

(2.10a)\begin{gather} \varPhi^{(11,0)} = \dfrac{gA_I(X-c_{g0}T)}{\omega_0}\dfrac{\cosh k_0(z+h_d)}{\cosh k_0 h_d} \sin\psi_0, \end{gather}

(2.10b)\begin{gather}\varPhi_I^{(11,1)}={-} \dfrac{g\partial_X A_I(X-c_{g0}T)}{\omega_0}\dfrac{(z+h_d)\sinh k_0(z+h_d)}{\cosh k_0 h_d} \cos\psi_0, \end{gather}

(2.10c)\begin{gather}\varPhi_I^{(22,0)} = \dfrac{3\omega_0A^2_I(X-c_{g0}T)}{8}\dfrac{\cosh 2k_0(z+h_d)}{\sinh^4 k_0 h_d} \sin 2\psi_0, \end{gather}

(2.10d)\begin{gather}\varPhi_I^{(20,1)} = \int\limits^{\kappa_m}_{-\kappa_m} \dfrac{\mathrm{i} g\omega_0\kappa\hat{E}(\kappa)B(k_0)}{g\kappa \tanh \delta\kappa h_d-\delta c^2_{g0}\kappa^2}\dfrac{\cosh \delta \kappa(z+h_d)}{\cosh \delta \kappa h_d} \exp({\mathrm{i} \kappa (X-c_{g0}T)})\,\mathrm{d} \kappa, \end{gather}

(2.10e)\begin{gather}\hat{E}(\kappa) = \dfrac{1}{4{\rm \pi}}\int^{\infty}_{-\infty} A_I^2(X-c_{g0}T)\, \mathrm{e}^{-\mathrm{i} \kappa X} \,\mathrm{d} X, \end{gather}

(2.10f)\begin{gather}B(k_0) = \dfrac{1}{\tanh k_0h_d} \left[\dfrac{c_{g0}}{2c_0}(1-\tanh^2k_0h_d)+1\right], \end{gather}

where $\kappa _m$ ($0<\kappa _m\ll k_0$) is the maximum wavenumber of the packet resulting from the assumption of narrow bandwidth, $c_0$ is the phase velocity and $c_{g0}$ is the group velocity of the wavepacket on the deeper side.

2.4. Overall structure of the solutions and underlying physics

Before constructing explicit solutions to the problem of interest, we first explain the key components of these solutions and the underlying physics. The solutions can be described as functions of the parameters of an incident wavepacket, as detailed in §§ 2.5 and 2.6. Taking the velocity potential as an example, a flow diagram of the solution associated with an incoming wavepacket is shown in figure 2, and a summary of the expressions for the velocity potential is presented in table 1 in Appendix D. In figure 2, the velocity potential is organised according to the order of product of wave steepness and bandwidth, as explained below. Naturally, we limit the discussion to those cases in which the incident wavepacket propagating over a step ‘feels’ the abrupt depth change. That is, the water depth compared to the carrier wavelength of an incoming wavepacket is ${O}(1)$ on at least one side of the step if not both.

Figure 2. Flow diagram of the perturbation theory solutions for the velocity potential of a narrow-banded wavepacket propagating over a step. The terms are organised according to the order of the product of wave steepness and bandwidth. From the top to the bottom row, the figure shows the incident, first-order, second-order superharmonic and second-order subharmonic or mean wavepackets. The subscripts $I$, $T$ and $E$ denote the incoming, transmitted and evanescent wavepackets, with $d$ and $s$ used to label the evanescent wavepackets on the deeper and shallower sides, respectively. The subscripts $b$ and $f$ denote bound and free waves at second order in wave steepness, respectively. A summary of the expressions for the velocity potential is given in Appendix D.

Table 1. Solutions for different velocity potentials at ${O}(\epsilon ^i\delta ^{\,j})$ according to the overall structure presented in § 2.4.

At first order in wave steepness, specifically ${O}(\epsilon \delta ^0)$, an incident wavepacket responds to an abrupt depth change by being reflected ($\varPhi ^{(11,0)}_R$) and transmitted ($\varPhi ^{(11,0)}_T$), complemented by the generation of evanescent waves ($\varPhi ^{(11,0)}_{Ed}$ on the deeper side and $\varPhi ^{(11,0)}_{Es}$ on the shallower side) near the step (cf. Massel 1983).

The mechanism that gives rise to waves at second order, namely ${O}(\epsilon ^2)$, can be divided into two parts. The first is the forcing of bound waves by combinations of linear waves that also arises in the absence of a step (cf. (2.10f)) and is well established (Massel 1983; Mei et al. 1989; Calvert et al. 2019). The second comprises the release of bound waves into free waves owing to the presence of the step. Forcing by combinations of linear waves leads to bound waves (denoted with the subscript $b$ in figure 2) that can only propagate together with the linear wavepacket. In contrast, free waves satisfy the linear dispersion relation and, hence, propagate independently. The bound waves include superharmonic bound waves (${O}(\epsilon ^2\delta ^0)$), which are proportional to $\exp (2\mathrm {i}\psi )$, and subharmonic bound waves (${O}(\epsilon ^2\delta ^1)$), which are independent of the rapidly varying phase $\psi _0$. Upon travelling over the step, these bound waves may be transmitted or reflected, staying bound, or be released into freely propagating wavepackets. In addition, new evanescent waves will be generated. The freely propagating wavepackets overlap with the linear wavepackets near the step, but will separate after a certain length of propagation owing to their different speeds. The distance over which separation occurs depends on the difference in group speeds and packet length.

2.5. First-order solutions (up to ${O}(\epsilon \delta ^1)$)

In this section, we extend the monochromatic-wave solutions presented in Massel (1983) to allow for a wavepacket that varies slowly in both space and time. Following Massel (1983), $\varPhi ^{(1)}$ is expressed as

(2.11a)\begin{align} \varPhi^{(1)} &= {\varPhi^{(11,0)}_I+\varPhi^{(11,0)}_R+\sum_{n=1}^{\infty}\varPhi^{(11,0)}_{Ed,n}} +\delta \varPhi^{(11,1)}+{O}(\epsilon\delta^2)\quad\textrm{for}\ x<0, \end{align}

(2.11b)\begin{align}\varPhi^{(1)} &= {\varPhi^{(11,0)}_T+\sum_{m=1}^{\infty}\varPhi^{(11,0)}_{Es,m}} + {\varPhi^{(11,1)}_T+\delta\varPhi^{(11,1)}}+{O}(\epsilon\delta^2)\quad\textrm{for}\ x>0, \end{align}

in which the subscripts $I$, $R$ and $T$ denote the (propagating) incoming, reflected and transmitted wavepackets, respectively. The subscripts $Ed,n$ and $Es,m$ denote the evanescent waves on the deeper and shallower sides, respectively. As for the case without a step, one can easily show that $\varPhi ^{(11,1)}$ does not contribute to the second-order solutions at ${O}(\epsilon ^2\delta )$, but only to those at higher orders in bandwidth (see § 2.6). The details of the derivation of $\varPhi ^{(11,1)}$ are nevertheless included in Appendix A for completeness.

The linearised boundary value problem (2.1)–(2.5) yields

(2.12a)\begin{gather} \varPhi^{(11,0)}_R =\dfrac{g}{\omega_0}\dfrac{\cosh k_0(z+h_d)}{\cosh k_0h_d}|R_0|A_R(X,T)\sin{({-}k_0x_0-\omega_0 t +\mu_0+\mu_R)} , \end{gather}

(2.12b)\begin{gather}\varPhi^{(11,0)}_T =\dfrac{g}{\omega_0}\dfrac{\cosh k_{0s}(z+h_s)}{\cosh k_{0s}h_s}|T_0|A_{T}(X,T)\sin{(k_{0s}x_0-\omega_0 t + \mu_0+ \mu_T)} , \end{gather}

(2.12c)\begin{gather}\varPhi^{(11,0)}_{Ed,n} = \mathcal{I} \left(\dfrac{g}{\omega_0}\dfrac{\cosh k_n(z+h_d)}{\cosh k_nh_d}R_n A_{Ed,n}(X,T) \exp({-\mathrm{i} k_nx_0-\mathrm{i}\omega_0 t+\mu_0})\right), \end{gather}

(2.12d)\begin{gather}\varPhi^{(11,0)}_{Es,m} = \mathcal{I} \left(\dfrac{g}{\omega_0}\dfrac{\cosh k_m(z+h_s)}{\cosh k_mh_s}T_m A_{Es,m}(X,T) \exp({\mathrm{i} k_mx_0-\mathrm{i}\omega_0 t+\mu_0})\right), \end{gather}

where the reflection coefficient $R_0=|R_0|\exp (\mathrm {i} \mu _R)$ and the transmission coefficient $T_0=|T_0|\exp (\mathrm {i} \mu _T)$ as well as the evanescent wave coefficients $R_n$ and $T_m$ are complex, and $\mathcal {I}$ denotes the imaginary component. The coefficients $R_0$, $R_n$, $T_0$ and $T_m$ of the free waves at ${O}(\epsilon \delta ^0)$ are solved for numerically based on the boundary conditions at the step described by (2.4) and (2.5):

(2.13a)\begin{align}\dfrac{\cosh k_0(z+h_d)}{\cosh k_0h_d} +\sum_{n=0}^{N} R_n \dfrac{\cosh k_n(z+h_d)}{\cosh k_nh_d} &= \sum_{m=0}^{M} T_m \dfrac{\cosh k_m(z+h_s)}{\cosh k_mh_s} \nonumber\\ &\qquad \textrm{for}\ -h_s < z < 0,\end{align}

(2.13b)\begin{align}\mathrm{i} k_0 \dfrac{\cosh k_0(z+h_d)}{\cosh k_0h_d} -\mathrm{i} k_n \sum_{n=0}^{N} R_n\dfrac{\cosh k_n(z+h_d)}{\cosh k_nh_d} &= \mathrm{i} k_m\sum_{j=0}^{M} T_m \dfrac{\cosh k_m(z+h_s)}{\cosh k_mh_s}\nonumber\\ &\qquad \textrm{for}\ -h_s < z < 0, \end{align}

(2.13c)\begin{align}\mathrm{i} k_0 \dfrac{\cosh k_0(z+h_d)}{\cosh k_0h_d} -\mathrm{i} k_n \sum_{n=0}^{N} R_n\dfrac{\cosh k_n(z+h_d)}{ \cosh k_nh_d} &= 0\quad \textrm{for}\ -h_d < z < {-}h_s, \end{align}

where $N$ and $M$ denote the finite number of evanescent modes used on the deeper and shallower sides, respectively. We show in Appendix B how $R_n$ and $T_m$ are numerically solved for using the orthogonality properties of the functions $\cosh k_i(z+h_d)$ and $\cosh k_i(z+h_s)$.

Departing from the analysis of Massel (1983), the packets are now allowed to vary slowly in time and space. Detailed derivations are presented in Appendix C. After taking into account the boundary conditions at the step, their dependence on time and space can be expressed as

(2.14a)\begin{gather} A_R(X,T) = A_I({-}X-c_{g0}T), \quad A_{Ed,n}(X,T) = A_I({c_{g0}}X/{c_{gn}}- c_{g0}T), \end{gather}

(2.14b)\begin{gather}A_T(X,T) = A_I({c_{g0}}X/{c_{g0s}}- c_{g0}T),\quad A_{Es,m}(X,T) = A_I({c_{g0}}X/{c_{gm}}- c_{g0}T), \end{gather}

where we note that these packets are continuous at $x=0$, the packets of the reflected and transmitted packets travel at the group speed determined by the local depth and we have used analytic continuation for the spatial dependence of the evanescent wavepackets. The following relations hold for the wavenumbers and group velocities, respectively:

(2.15)\begin{gather} \omega^2_0=gk_i\tanh k_ih_d = gk_j \tanh k_j h_s, \end{gather}

(2.16a,b)\begin{gather} c_{gi} =\left. \dfrac{\text{d}\omega}{\text{d} k}\right|_{k=k_{i}} ,\quad c_{gj} =\left. \dfrac{\text{d}\omega}{\text{d} k}\right|_{k=k_{j}}, \end{gather}

where $i = {0}$ or ${n}$, $j = {0s}$ or ${m}$, $k_0$ and $k_{0s}$ are real wavenumbers, and the rest of the wavenumbers and corresponding group speeds are imaginary. The imaginary wavenumbers correspond to evanescent waves that vanish with horizontal distance away from the step, as $\exp (-|\mathrm {i} k_n x|)$ or $\exp (-|\mathrm {i} k_m x|)$.

2.6. Second-order solutions $({O}(\epsilon ^2))$

The free-surface boundary conditions can be combined into one, which gives at ${O}(\epsilon ^2)$ (e.g. Longuet-Higgins & Stewart 1964; McAllister et al. 2018)

(2.17)\begin{align} \dfrac{\partial^2 \varPhi^{(2)}}{\partial t^2} + g \dfrac{\partial \varPhi^{(2)}}{\partial z} &={-} \frac{\partial}{\partial t} \left[ \frac{1}{2}\left(\dfrac{\partial \varPhi^{(1)}}{\partial x} \right)^2+ \frac{1}{2}\left(\dfrac{\partial \varPhi^{(1)}}{\partial z} \right)^2 + \zeta^{(1)}\frac{\partial^2 \varPhi^{(1)}}{\partial z\partial t} \right] \nonumber\\ &\quad + g \dfrac{\partial}{\partial x} \left( \frac{\partial \varPhi^{(1)}}{\partial x}\zeta^{(1)}\right) \quad \text{at}\ z=0,~ \end{align}

with a corresponding diagnostic equation for the surface elevation at ${O}(\epsilon ^2)$

(2.18)\begin{equation} \zeta^{(2)} ={-}\dfrac{1}{g}\left(\dfrac{\partial \varPhi^{(2)}}{\partial t} + \dfrac{1}{2} \left(\boldsymbol{\nabla}\varPhi^{(1)} \right)^2+\dfrac{\partial^2 \varPhi^{(1)}}{\partial t\partial z}\zeta^{(1)}\right), \end{equation}

in which the second term on the right-hand side becomes zero in deep water.

Substituting the linear solutions (2.11) into (2.17) and collecting terms at ${O}(\epsilon ^2)$ yields

(2.19)\begin{equation} \dfrac{\partial^2 \varPhi^{(2)}}{\partial t^2} + g \dfrac{\partial \varPhi^{(2)}}{\partial z} = \delta^0 F_{\epsilon^2\delta^0} +\delta F_{{\epsilon^2\delta^1}} +{O}(\epsilon^2\delta^2) \end{equation}

with

(2.20)\begin{align} F_{\epsilon^2\delta^0} &={-} \frac{\partial}{\partial {t_0}}\left[\frac{1}{2}\left(\dfrac{\partial \varPhi^{(11,0)}}{\partial {x_0}} \right)^2+ \frac{1}{2}\left(\dfrac{\partial \varPhi^{(11,0)}}{\partial z} \right)^2 + \zeta^{(11,0)}\frac{\partial^2 \varPhi^{(11,0)}}{\partial z\partial {t_0}}\right] \nonumber\\ &\quad + g\dfrac{\partial}{\partial {x_0}} \left( \frac{\partial \varPhi^{(11,0)}}{\partial {x_0}}\zeta^{(11,0)}\right), \end{align}

(2.21)\begin{align} F_{\epsilon^2\delta^1} &={-} \frac{\partial}{\partial {T}} \left[ \frac{1}{2}\left(\dfrac{\partial \varPhi^{(11,0)}}{\partial {x_0}} \right)^2+ \frac{1}{2}\left(\dfrac{\partial \varPhi^{(11,0)}}{\partial z} \right)^2 + \zeta^{(11,0)}\frac{\partial^2 \varPhi^{(11,0)}}{\partial z\partial {t_0}}\right]\nonumber\\ &\quad - \frac{\partial}{\partial {t_0}} \left[ \dfrac{\partial \varPhi^{(11,0)}}{\partial {x_0}} \left( \dfrac{\partial \varPhi^{(11,0)}}{\partial {X}}+ \dfrac{\partial \varPhi^{(11,1)}}{\partial {x_0}}\right) + \dfrac{\partial \varPhi^{(11,0)}}{\partial {z}} \dfrac{\partial \varPhi^{(11,1)}}{\partial {z}}\right]\nonumber\\ &\quad +\frac{\partial}{\partial {t_0}} \left( \zeta^{(11,1)}\frac{\partial^2 \varPhi^{(11,0)}}{\partial z\partial {t_0}} + \zeta^{(11,0)}\frac{\partial^2 \varPhi^{(11,1)}}{\partial z\partial {t_0}} \right) \nonumber\\ &\quad + g \dfrac{\partial}{\partial {X}} \left(\frac{\partial \varPhi^{(11,0)}}{\partial {x_0}}\zeta^{(11,0)}\right) + g \dfrac{\partial}{\partial {x_0}} \left(\zeta^{(11,1)}\frac{\partial \varPhi^{(11,0)}}{\partial {x_0}} + \zeta^{(11,0)}\frac{\partial \varPhi^{(11,1)}}{\partial {x_0}} \right), \end{align}

in which $F_{\epsilon ^2\delta ^0}$ and $F_{\epsilon ^2\delta ^1}$ can be further decomposed based on wave harmonics. A similar equation can be obtained for $\zeta ^{(2)}$ (not shown).

After identifying the harmonics of the forcing terms $F_{\epsilon ^2\delta ^0}$ and $F_{\epsilon ^2\delta ^1}$, the second-order solutions of $\varPhi ^{(2)}$ and $\zeta ^{(2)}$ can be separated into two parts with different harmonics:

(2.22a)\begin{gather} \varPhi^{(2)} = \varPhi^{(22,0)} + \delta \varPhi^{(20,1)}, \end{gather}

(2.22b)\begin{gather}\zeta^{(2)}= \zeta^{(22,0)}+\delta \zeta^{(20,1)}, \end{gather}

in which $\varPhi ^{(22,0)}$ and $\zeta ^{(22,0)}$ are the superharmonic terms proportional to $\exp {(-2\mathrm {i}\omega _0t_0)}$ and $\varPhi ^{(20,1)}$ and $\zeta ^{(20,1)}$ are the subharmonic (or mean) terms that are independent of fast time $t_0$.

2.6.1. Superharmonic packets at ${O}(\epsilon ^2\delta ^0)$

Similar to Massel (1983), we seek solutions for the superharmonic packets at ${O}(\epsilon ^2\delta ^0)$ of the form

(2.23)\begin{gather} \varPhi^{(22,0)} = \varPhi^{(22,0)}_I + \varPhi^{(22,0)}_{R,b} + \varPhi^{(22,0)}_{R,f} + \mathcal{I}\left(\sum_{n=1}\varPhi^{(22,0)}_{Ed,n}\right)\quad\text{for}\ x<0, \end{gather}

(2.24)\begin{gather}\varPhi^{(22,0)} = \varPhi^{(22,0)}_{T,b} + \varPhi^{(22,0)}_{T,f} +\mathcal{I}\left( \sum_{m=1}\varPhi^{(22,0)}_{Es,m}\right)\quad\text{for}\ x>0, \end{gather}

in which subscripts $b$ and $f$ denote superharmonic bound and free waves, respectively. In order to obtain tractable solutions, we ignore forcing by product linear evanescent waves, which are typically small. We justify this assumption ex post by comparing with fully nonlinear numerical simulations. We note that inclusion of forcing by evanescent terms can lead to convergence problems of second-order solutions (Monsalve Gutiérrez 2017). After considerable manipulation, we obtain

(2.25a)\begin{gather} \varPhi^{(22,0)}_{R,b} = \dfrac{3\omega_0}{8}|R_0|^2A^2_R\dfrac{\cosh 2k_0(z+h_d)}{\sinh^4 k_0 h_d} \sin({-}2k_0x-2\omega_0 t + 2\mu_0 +2\mu_R), \end{gather}

(2.25b)\begin{gather}\varPhi^{(22,0)}_{R,f} = \omega_0|R_{20}|A^2_I\left(-\dfrac{c_{g0}}{c_{g20}}X-c_{g0}T\right)\dfrac{\cosh k_{20}(z+h_d)}{\cosh k_{20} h_d} \sin({-}k_{20}x-2\omega_0 t + 2\mu_0+ \mu_{20}), \end{gather}

(2.25c)\begin{gather}\varPhi^{(22,0)}_{Ed,n} = \omega_0 R_{2n} A^2_I\left(-\dfrac{c_{g0}}{c_{g2n}}X-c_{g0}T\right)\dfrac{\cosh k_{2n}(z+h_d)}{\cosh k_{2n} h_d} \exp(-\mathrm{i} (k_{20}x+2\omega_0 t) + 2\mu_0), \end{gather}

(2.25d)\begin{gather}\varPhi^{(22,0)}_{T,b} = \dfrac{3\omega_0}{8}|T_0|^2A^2_I\left(\dfrac{c_{g0}}{c_{g0s}}X-c_{g0}T\right)\dfrac{\cosh 2k_{0s}(z+h_s)}{\sinh^4 k_{0s} h_s} \sin(2k_{0s}x-2\omega_{0s} t + 2\mu_0 + 2\mu_T), \end{gather}

(2.25e)\begin{gather}\varPhi^{(22,0)}_{T,f} = \omega_0 |T_{20}|A^2_I\left(\dfrac{c_{g0}}{c_{g20s}}X-c_{g0}T\right)\dfrac{\cosh k_{20s}(z+h_s)}{\cosh k_{20s} h_s}\sin(k_{20s}x-2\omega_0 t + 2\mu_0+ \mu_{20s}), \end{gather}

(2.25f)\begin{gather}\varPhi^{(22,0)}_{Es,m}= \omega_0 T_{2m} A^2_I\left(\dfrac{c_{g0}}{c_{g2m}}X-c_{g0}T\right) \dfrac{\cosh k_{2m}(z+h_s)}{\cosh k_{2m} h_s} \exp(-\mathrm{i} (k_{2m}x-2\omega_0 t + 2\mu_0)) , \end{gather}

(2.25g)\begin{gather}4\omega_0^2 = g k_{2i} \tanh (k_{2i}h_d)\quad \textrm{for}\ i=0, 1, 2, 3,\ldots, \end{gather}

(2.25h)\begin{gather}4\omega_0^2 = g k_{2j} \tanh (k_{2j} h_s)\quad \textrm{for}\ j=0s, 1, 2, 3,\ldots, \end{gather}

where the last two equations denote the dispersion relationships associated with the frequency $2\omega _0$ for two different depths $h_d$ and $h_s$. The wavenumbers $k_{20}$ and $k_{20s}$ are real, and the other superharmonic wavenumbers are imaginary and correspond to evanescent waves. The reflection ($R_{2n}$) and transmission ($T_{2m}$) coefficients are solved for numerically from the boundary conditions at the step (2.4) and (2.5) at this particular order. The group velocities and the phases are defined as

(2.26a)\begin{gather} c_{g2i} = \dfrac{\omega_0}{k_{2i}}\left(1+\dfrac{2k_{2i}h_j}{\sinh(2k_{2i}h_j)}\right), \end{gather}

(2.26b)\begin{gather}\psi_{20} = \text{arg}(R_{20}),\quad \psi_{20s} = \text{arg}(T_{20}). \end{gather}

in which $k_{2i}h_j = k_{2i}h_d$ ($i=0,1,2,\ldots$) or $k_{2i}h_j = k_{2i}h_s$ ($i=0s,1,2,\ldots$).

2.6.2. Subharmonic packets at ${O}(\epsilon ^2\delta ^1)$

In this section, we present second-order subharmonic solutions, which were not included by Massel (1983). Averaging in time over the fast scales, we find that $F_{\epsilon ^2\delta ^0}=0$. The solutions at ${O}(\epsilon \delta )$ do not contribute to leading order, and we obtain the following for the subharmonic forcing at second order:

(2.27a)\begin{gather} \overline{F_{\epsilon^2\delta^0}+\delta F_{\epsilon^2\delta^1}}={-} \delta \dfrac{g^2}{\omega^2_0}\frac{\partial}{\partial {T}}\left[\frac{k_0^2}{4}(1- \tanh^2k_0h_d)(A_I^2+A_R^2)\right] \nonumber\\ \hspace{5.3pc}+\,\delta \dfrac{g^2}{2\omega_0} \frac{\partial}{\partial {X}} [k_0(A_I^2-A_R^2)] \quad \text{for}\ X<0, \end{gather}

(2.27b)\begin{gather} \overline{F_{\epsilon^2\delta^0}+\delta F_{\epsilon^2\delta^1}}={-} \delta \dfrac{g\omega_0}{2c_{g0s}\tanh k_{0s}h_s} \left[\dfrac{c_{g0s}}{2c_{0s}}(1-\tanh^2k_{0s}h_s)+1\right]\partial_T(A_T^2)\quad \text{for}\ X>0, \end{gather}

where $c_{0s}=\omega _0/k_{0s}$ denotes the phase velocity of the linear carrier wave on the shallower side. In order to maintain tractable solutions, we ignore forcing by linear evanescent waves. As for the second-order superharmonic solutions, we justify this assumption ex post by comparing with fully nonlinear numerical simulations.

The forcing in (2.27), together with the Laplace equation (2.1) and the bottom boundary condition (2.3), leads to the following subharmonic bound waves:

(2.28a)\begin{align} \varPhi^{(20,1)}_b(X,T,z) &= \varPhi_I^{(20,1)} + B_R {\int\limits_{-\varOmega_m}^{\varOmega_m}} \dfrac{\mathrm{i}\kappa g\omega_0\hat{E}_R(\varOmega)}{g\kappa_0\tanh \delta \kappa_0 h_d-\delta \varOmega^2} \dfrac{\cosh \delta \kappa_0(z+h_d)}{\cosh \delta \kappa_0 h_d}\, \mathrm{e}^{-\mathrm{i} \varOmega T}\,\mathrm{d} \varOmega \nonumber\\ &\qquad\text{for}\ X<0,\end{align}

(2.28b)\begin{align} \varPhi^{(20,1)}_b(X,T,z) &= B_T{\int\limits_{-\varOmega_m}^{\varOmega_m} \dfrac{\mathrm{i} {\kappa_{0s}} g \omega_0\hat{E}_T({\varOmega})}{g{\kappa_{0s}}\tanh\delta{\kappa_{0s}} h_s-\delta \varOmega^2} } \dfrac{\cosh \delta \kappa_{0s}(z+h_s)}{\cosh \delta \kappa_{0s} h_s}\,\mathrm{e}^{-\mathrm{i} \varOmega T}\,\mathrm{d} \varOmega \nonumber\\ &\qquad \text{for}\ X>0, \end{align}

where $\kappa _0 = -\varOmega /c_{g0}$, $\kappa _{0s} = -\varOmega /c_{g0s}$, $\varOmega _m$ ($0<\delta \varOmega _m\ll \omega _0$) is the maximum frequency of the packet resulting from the assumption of narrow bandwidth and

(2.29a)\begin{gather} \hat{E}_R(\varOmega) = \dfrac{1}{4{\rm \pi}} \int\limits^{\infty}_{-\infty} R_0^2A_R^2\, \mathrm{e}^{\mathrm{i}\varOmega T}\,\mathrm{d} T, \quad \hat{E}_T{(\varOmega)} = \dfrac{1}{4{\rm \pi}} \int\limits^{\infty}_{-\infty} T_0^2A_T^2\,\mathrm{e}^{\mathrm{i}\varOmega T}\,\mathrm{d} T, \end{gather}

(2.29b)\begin{gather}B_R ={-} \dfrac{1}{\tanh k_0h_d} \left[\dfrac{c_{g0}}{2c_0}(1-\tanh^2k_0h_d)-1\right], \end{gather}

(2.29c)\begin{gather}B_T = \dfrac{1}{\tanh k_{0s}h_s} \left[\dfrac{c_{g0s}}{2c_{0s}}(1-\tanh^2k_{0s}h_s)+1\right]. \end{gather}

The bound subharmonic waves in (2.28) correspond to those of the incoming, reflected and transmitted separately. Together, these bound waves do not satisfy the boundary conditions at the step, where additional free waves are generated. To avoid prohibitively cumbersome solutions, we make the additional assumption that the subharmonic packet is long relative to the water depth, so that the mean flow is shallow (see Calvert et al. 2019). This assumption covers most practical applications in coastal waters.

2.6.3. The long-wave approximation for subharmonic packets ($1/(kh\delta ) \gg 1$)

In the limit $1/(kh\delta ) \gg 1$ for both $h_d$ and $h_s$, but for $k_0 h_d={O}(1)$ and $k_{0s} h_s={O}(1)$, the bound subharmonic behaviour can be described in terms of horizontal velocities

(2.30a–c)\begin{equation} u^{(20,1)}_{I} = \dfrac{g k_0 B_d}{c_{g0}}A_I^2,\quad u^{(20,1)}_{R,b} ={-} \dfrac{g k_0|R_0|^2 B_d}{c_{g0}}A_R^2,\quad u^{(20,1)}_{T,b} =\dfrac{g k_{0s}|T_0|^2 B_s }{c_{g0s}} A_T^2, \end{equation}

and mean set-downs of the surface elevation

(2.31a–c)\begin{equation} \zeta^{(20,1)}_{I} =k_0 B_d A_I^2, \quad \zeta^{(20,1)}_{R,b} =k_0B_d |R_0|^2 A_R^2,\quad \zeta^{(20,1)}_{T,b} =k_{0s} B_s |T_0|^2 A_T^2, \end{equation}

where the non-dimensional coefficients $B_d$ and $B_s$ are given by

(2.32a)\begin{gather} B_d ={-}\dfrac{1}{4 (g h_d- {c}_{g0}^2 )} \left(\dfrac{2gh_d - {c}_{g0} ^2 }{2 \sinh 2k_0h_d} +\dfrac{2g c_{g0} }{\omega_0} \right) , \end{gather}

(2.32b)\begin{gather}B_s ={-} \dfrac{1}{4 (g h_s- {c}_{g0s}^2 )} \left(\dfrac{2gh_s - {c}_{g0s} ^2 }{2\sinh 2k_{0s}h_s} +\dfrac{2g c_{g0s} }{\omega_0} \right). \end{gather}

When the carrier waves are additionally assumed to travel in deep water (i.e. $k_0h_d \gg 1$ and $k_{0s}h_s\gg 1$) then the first term in the brackets of both (2.32a) and (2.32b) vanishes. For completeness, we note that, owing to limit $h/\sigma \to 0$, the order of the solutions in $\delta$ has increased by one, although we do not update our notation to reflect this.

In accordance with the long-wave approximation for the subharmonic packets, freely travelling subharmonic packets generated at the step propagate at the shallow-water velocity, i.e. $\sqrt {gh_d}$ on the deeper side and $\sqrt {gh_s}$ on the shallower side. Assuming such free subharmonic packets can propagate in both directions, we seek solutions of the form

(2.33a)\begin{gather} \zeta^{(20,1)}_{R,f} = B_R^f k_0 A^2_I\left(-\frac{c_{g0}}{\sqrt{gh_d}}X-c_{g0}T\right), \quad \zeta^{(20,1)}_{T,f} = B_T^f k_{0s}A^2_I\left(\frac{c_{g0}}{\sqrt{gh_s}}X-c_{g0}T\right), \end{gather}

(2.33b)\begin{gather} u^{(20,1)}_{R,f} ={-} \sqrt{\frac{g}{h_d}}B_R^f k_{0} A^2_I\left(-\frac{c_{g0}}{\sqrt{gh_d}}X-c_{g0}T\right), \quad u^{(20,1)}_{T,f} = \sqrt{\frac{g}{h_s}}B_T^f k_{0s} A^2_I\left(\frac{c_{g0}}{\sqrt{gh_s}}X-c_{g0}T\right). \end{gather}

The relationship between $u$ and $\zeta$ is set by $\partial _t u = -g \partial _x \zeta$ (cf. (4.1.3) in Mei et al. (1989)). The coefficients $B_R^{f}$ and $B_T^{f}$ must be found from the matching conditions at the step. For a shallow flow (see e.g. Mei et al. (1989) for details), these become (i) continuity of the volume flux across the step and (ii) continuity of the free surface across the step:

(2.34a,b)\begin{equation} \lim_{X\to0^-}u^{(20,1)} h_d = \lim_{X\to0^+}u^{(20,1)} h_s,\quad\lim_{X\to0^-}\zeta^{(20,1)} =\lim_{X\to0^+}\zeta^{(20,1)}, \end{equation}

where we note that $u^{(1)}(z=0)\eta ^{(1)}$ from depth integration of the linear velocity truncated at second order is not included in (2.34a), as this is already continuous across the step. Hence, we obtain

(2.35a)\begin{gather} B_{T}^f = \dfrac{gk_{0s}h_sB_s|T_0|^2/c_{g0s} + gk_{0}h_dB_d(1-|R_0|^2)/c_{g0}+ \sqrt{gh_d} (k_{0s}B_s|T_0|^2 - k_{0}B_d(|R_0|^2+1))} {(\sqrt{gh_d}+\sqrt{gh_s})k_{0s}}, \end{gather}

(2.35b)\begin{gather}B_{R}^f = \dfrac{gk_{0s}h_sB_s|T_0|^2/c_{g0s} - gk_{0}h_dB_d(1-|R_0|^2)/c_{g0} - \sqrt{gh_s} (k_{0s}B_s|T_0|^2 - k_{0}B_d(|R_0|^2+1))} {(-\sqrt{gh_d}-\sqrt{gh_s})k_{0}}. \end{gather}

We note that the relations $\text {sign}(B_T^f) = - \text {sign}(B_s)$ and $\text {sign}(B_R^f) = - \text {sign}(B_s)$ hold, and both free waves are thus positive, taking the form of set-ups, as the sign of the bound set-down is always negative. The coefficients $B_{T}^f$ and $B_{T}^f$ only depend on two non-dimensional parameters: $k_0h_d$ and $k_0h_s$. We further explore these solutions and the underlying physics in the next section.

3. Results

In order to examine the predictions of the theoretical model in § 2, we consider an incoming Gaussian wavepacket on the deeper side defined as follows:

(3.1)\begin{equation} \zeta^{(11,0)}_I = A_0 \exp\left(-\dfrac{(x-x_\text{f} - c_{g0}(t-t_{f}))^2}{2\sigma^2_x}\right )\cos(k_0x - \omega_0t), \end{equation}

in which $k_0$ and $\omega _0$ are the carrier wavenumber and angular frequency, respectively, $c_{g0} = \omega _0 /(2k_0) (1+2k_0h_d/\sinh (2k_0h_d))$ is the group velocity on the deeper side, $\sigma _x$ is the characteristic length of the packet, and $x_{f}$ denotes the location where the linear irregular waves focus at time $t_{f}$. We set the wave steepness $\epsilon = k_0A_0=0.03$ and the bandwidth parameter $\delta = 1/(k_0\sigma _x)=0.06$, so that both remain much smaller than $1$ in accordance with the assumptions presented in § 2.

We examine three distinct stages of evolution: stage I when the packet is sufficiently far ahead of the step on the deeper side, stage II when the packet ‘feels’ the step and transient processes in the vicinity of the step take place and stage III when the packet has left the step behind.

3.1. Generation of free packets: stage I versus stage III

Figure 3 shows the theoretically predicted free-surface elevation before (stage I) and after (stage III) passing the step. Before the step ( figure 3a–d), the main (linear) packet is associated with an in-phase superharmonic bound wavepacket and a subharmonic bound set-down (cf. (2.10f)), as is well known (e.g. Mei et al. 1989; Calvert et al. 2019). After the step ( figure 3e–h), both the superharmonic bound wavepacket and the subharmonic bound set-down have increased in magnitude. Also present are two additional superharmonic wavepackets and two additional subharmonic components, only one of which is visible in figure 3.

Figure 3. Theoretically predicted interaction with a step decomposed by order and phase, showing stage I before the step (a–d) and stage III after the step (e–h). In the figure, $\epsilon = 0.03$, $\delta = 0.06$, $k_0h_d = 1.1$, $k_{0s}h_s = 0.70$, $h_d = 0.75$ m, $h_s/h_d = 0.53$, $T_{{p}} = 1.9$ s is the carrier wave period and $\lambda _0 = 4.4$ m. The step is located at $x = 0$. Panels (a–d) correspond to $t = 30 T_p$ before the main wavepacket reaches the step and (e–f) to $t = 70 T_p$ after it has passed the step.

The response to the step is most clearly illustrated in figure 4. Focusing on figure 4(a) first, the bound superharmonic wavepacket on the deeper side is split into three wavepackets after experiencing the depth transition, one of which stays bound and travels with the main packet at $c_{g0s}$. A first additional superharmonic free wavepacket propagates in the same direction as the main packet, but slower at $c_{g20}$ ($c_{g20s} < c_{g0s}$). A second additional superharmonic free wavepacket propagates and is reflected and travels in the opposite direction at an absolute speed of $c_{g20}$ ($c_{g20} < c_{g0}$).

Figure 4. Theoretically predicted spatio-temporal evolution of superharmonic (a) and subharmonic (b) wavepackets following interaction with a step at $x=0$. The parameters are the same as in figure 3. The straight red lines with arrowheads indicate the different group speeds and their propagation directions.

Analogous behaviour is observed in figure 4(b), except that the subharmonic free components are shallow-water waves and travel at higher speeds than the main (linear) packet. The subharmonic bound wave, manifest as a set-down of the free surface, becomes deeper on the shallower side. A free subharmonic set-up is released that propagates at the shallow-water speed $\sqrt {gh_s}$ in the direction of the main packet but faster. A free subharmonic set-down is reflected and travels in the opposite direction to the main packet at an absolute speed $\sqrt {gh_d}$ on the deeper side.

3.2. Amplitudes change and phases shift due to an abrupt depth transition

In the previous section, we examined a single combination of parameters. Although the four additional free second-order components will be generated for any combination of parameters, their amplitudes and phases depend on two dimensionless parameters, the relative depth on the deeper side $k_0h_d$ and the depth ratio $h_d/h_s$, in addition to the steepness squared. The linear reflection and transmission coefficients $R_0$ and $T_0$ are computed based on (2.13). The reflection and transmission coefficients for the superharmonic free wavepackets $R_{20}$ and $T_{20}$ are computed from the boundary conditions at the step in a similar fashion to (2.13). The reflection and transmission coefficients for the subharmonic free components $B_R^f$ and $B_T^f$ are given explicitly in (2.34a,b). The magnitudes and phases of all these coefficients are shown as functions of the two non-dimensional parameters in figure 5.

Figure 5. Contours of amplitudes and phases of the theoretically predicted reflection and transmission coefficients as functions of $k_0h_d$ and the depth ratio $h_s/h_d$.

Examining first the coefficients for the first harmonic shown in figure 5(a–d) (see also Massel 1983), the transmitted waves are amplified for $k_0h_d \lesssim 2.0$. The coefficient of reflection can reach a maximum of $\sim 30\,\%$ when the depth ratio decreases to $0.3$. Figure 5(b,d) shows that, relative to the incoming wavepacket, the transmitted linear waves generally have small phase shifts (${\lesssim }0.05{\rm \pi}$) and the reflected waves have a phase shift ${\lesssim }0.2{\rm \pi}$ when their amplitudes are ${\sim }10\,\%$–$30\,\%$ of the incoming wavepacket (comparing figures 5a and 5b).

Figure 5(e, f,i,j) shows the reflection and transmission coefficients of the free superharmonic waves. These are generally largest in magnitude for small $k_0 h_d$ and small depth ratios $h_s/h_d$ with the transmitted component considerably larger than the reflected component. Relative to the incoming wavepacket, the reflected superharmonic waves show small phase shifts, whereas the transmitted waves have a phase shift of between $-0.9{\rm \pi}$ and $-{\rm \pi}$. The latter is the cause of local transient maxima in crest elevation occurring in the vicinity of the step, as we examine in § 3.3.

The coefficients for the reflected and transmitted free and bound subharmonic components are calculated based on the long-wave approximation for these components presented in § 2.6.3 and shown in figure 5(g,h,k,l). The reflected free subharmonic components travel backwards on the deeper side in the form of a set-down, whereas the transmitted free subharmonic components travel forwards on the shallower side in the form of a set-up. We emphasise that the coefficients presented in figure 5 need to be used with care for small $k_0 h_s$, as a Stokes expansion is likely no longer valid for very shallow depths.

3.3. Behaviour near the abrupt depth transition: stage II

As noted in § 3.2, the transmitted superharmonic free wavepacket has a phase shift of approximately ${\rm \pi}$ relative to the transmitted main wavepacket (and its in-phase bound superharmonics). As a result of the smaller group velocity, the superharmonic and the transmitted main packet temporarily overlap just after the step before separating. These processes can be associated with two characteristic length scales: a beating length $L_b$ and an overlapping length $L_o$. Beating occurs when the free and bound superharmonic waves are in phase, namely at $x=L_b$ and

(3.2)\begin{equation} L_b(n)=\dfrac{-\arg(T_{20})+2\arg(T_{0})+2{\rm \pi} (n-1)}{k_{20s}-2k_{0s}}, \end{equation}

for any positive integer $n$ with the first beat corresponding to $n=1$, noting that $\arg (T_{0})\approx 0$ and $\arg (T_{2})\approx -{\rm \pi}$. Taking $4\sigma _{x,s}$ with $\sigma _{x,s}=\sigma _x c_{g0s}/c_{g0}$ as an estimate of the length of the group, the two groups will no longer overlap at $x=L_o$ and

(3.3)\begin{equation} L_o=\dfrac{4\sigma_{x,s}}{1-c_{g20s}/c_{g0s}}, \end{equation}

which denotes the distance between the peak of the main wavepacket and the step when the two groups just separate.

As the envelopes of the superharmonic bound and the superharmonic free waves travel at different group speeds and the lengths of the packets are limited, observation of $n$ beats requires $L_b(n)\leqslant L_o$. Hence, only the first (few) beat(s) will be observed. The length scales $L_o$ and $L_b(n)$ scaled by the carrier wavelength on the shallower side $\lambda _{0,s}$ are shown in figure 6 as a function of the dimensionless depth $k_{0s}h_s$. We can observe from figure 6(a) that the length for the first beat increases rapidly as the shallower water depth $k_{0s}h_s$ decreases for $k_{0s}h_s\lesssim 1.5$. At least one beat can be expected for $k_{0s}h_s>0.2$ and $\delta _{s}=1/(k_{0s}\sigma _{x,s})<0.1$, as shown in figure 6(b). As the group length increases (i.e. $\delta _{s}$ decreases), more beats can be expected. When we eventually approach to the limit $\delta _{s} \to 0$ (not shown here), denoting a uniform Stokes wave as studied and examined in Massel (1983), there are an infinite number of beats.

Figure 6. (a) The beating length ($L_b(n)$ for $n=1$) and (b) the overlapping and beating length ($L_o$ and $L_b(n)$ for $n\geqslant 1$) as a function of $k_{0s}h_s$. In the figure, $\lambda _{0,s} = 2{\rm \pi} /k_{0s}$ is the carrier wavelength of the main packet on the shallower side; $\delta _s = 1/(k_{0s}\sigma _{x,s})$ denotes the non-dimensional bandwidth on the shallower side.

4. Numerical validation

4.1. A fully nonlinear potential-flow numerical solver

In order to validate our solutions and justify our assumption that evanescent waves do not contribute meaningfully to behaviour at second order in steepness when waves travel over a step, we perform fully nonlinear potential-flow simulations. We employ a fully nonlinear potential-flow numerical solver that uses the boundary element method for the boundary value problem described by (2.1)–(2.5). The resulting numerical wave tank was first developed by Koo & Kim (2004) and has recently been used to examine a related problem by Zheng et al. (2020). Generation of waves in this numerical wave tank is based on linear theory (Havelock 1929), consistent with our experiments reported on in the companion paper (Li et al. 2021).

4.2. Comparison between theory and numerical simulations

Although we have tested a number of different cases, we show a comparison here between theory and numerical simulations for the same parameters as examined in the previous section, namely $\epsilon = 0.03$, $\delta = 0.06$, $k_0h_d = 1.1$, $k_{0s}h_s = 0.70$, $h_d = 0.75$ m, $h_s/h_d = 0.53$, $T_p = 1.9$ s and $\lambda _0 = 4.4$ m, which denotes the carrier wavelength on the deeper side. We choose a Gaussian wavepacket that focuses linearly at $t_{f} = 30T_p$ and $x_{{f}}=-5\lambda _0$ with the step at $x=0$. A computational domain of a length equal to $100\lambda _0$ was chosen in the numerical wave tank, where the distances between the wavemaker and the step and between the damping zone (the beach) are equal to $20\lambda _0$ and $80\lambda _0$, respectively. In order to compare the theoretical and numerical solutions explicitly for each order and phase (sub- and superharmonic), we filter the (narrow-banded) numerical solutions in the frequency domain.

Figure 7, which shows a comparison of the theoretically and numerically predicted surface elevations near the step decomposed by order and harmonic, reveals almost perfect agreement for the first-order and the second-order superharmonic and subharmonic surface elevation. The agreement includes phase and amplitude in the vicinity of the step. Figure 8 confirms this good agreement in Fourier space and figure 9 in physical space. The agreement becomes less perfect in the long-time tail of the wavepackets in figure 7, reflecting the non-dispersive approximation made in our theoretical solutions owing to truncation in bandwidth and nonlinearity (cf. Tayfun 1980, 1986; Mei et al. 1989; Trulsen & Dysthe 1996). The almost negligible difference between the theory and numerical simulations shown in figure 7(i) is due to the long-packet approximation we made in § 2.6.3.

Figure 7. Comparison of the theoretically and numerically predicted interaction with the step decomposed by order. The wavepacket focuses linearly at $t= 30 T_p$ and $x=-5\lambda _0$ with the step at $x=0$. The three different times ($t= 30 T_p$, $35T_p$, $40T_p$) corresponding to the three columns capture the passing of the step.

Figure 8. Comparison of the theoretically and numerically predicted free-surface elevation in Fourier wavenumber space ($k$) at different times (a–c) and in frequency space ($f$) at different locations (d–f). The parameters are the same as in figure 7.

Figure 9. Comparison of the theoretically (red dashed lines) and numerically (black solid lines) predicted spatio-temporal evolution of superharmonic (a) and subharmonic (b) wavepackets following interaction with a step at $x=0$ (vertical dashed line). The parameters are the same as in figure 7. The straight blue lines with arrowheads indicate the different group speeds and their propagation directions. The numerical results are represented by solid black lines and theory by dashed red lines.

5. Conclusions

This paper has examined the interaction of deterministic surface gravity wavepackets with an abrupt depth transition in the form of a step in intermediate water depth. To do so, we have developed second-order theory for narrow-banded wavepackets based on a Stokes and multiple-scales expansion, thereby extending the work of Massel (1983), which is only valid for monochromatic waves. To obtain tractable solutions from perturbation theory, we additionally assume that forcing of second-order terms due to local first-order evanescent waves, which are generated at first order owing to the presence of the step and vanish rapidly with distance away from it, can be ignored. We justify this assumption ex post by performing numerical simulations using a fully nonlinear potential-flow solver. Good agreement with our theoretical solutions is found.

As a wavepacket travels over a sudden depth transition, additional wavepackets are generated that propagate freely obeying the linear dispersion relation and arise at both first and second order in wave steepness in a Stokes expansion. As the superharmonic bound waves travel over the step, their magnitude changes, and two freely travelling superharmonic wavepackets are released. The two free packets consist of a generally in-phase reflected packet that travels in the opposite direction, and a generally out-of-phase transmitted packet propagating in the same direction as the main (linear) packet albeit at a lower speed. The same happens for the subharmonic components. In the subharmonic shallow-water limit, in which the packet is long relative to the water depth, we can find solutions for the subharmonic components in closed form. At the subharmonic level, the bound set-down generally increases in magnitude, and a free transmitted set-up travels ahead of the main packet with a free set-down being reflected.

In the region near the top of a depth transition, the resulting transient processes play a crucial role. In intermediate water depths, these processes are generally dominated by the superharmonic terms. Both the superharmonic bound waves and the freely travelling superharmonic waves appear immediately after the step on the shallower side. This causes beating, which is modulated by the envelopes of both packets and only exists near the step. Together, these effects cause a series of local peaks in surface elevation, which decline in magnitude with distance away from the step, declining more strongly for less narrow-banded and thus spatially shorter wavepackets. Such a decline is absent in the spatially periodic beating pattern predicted by Massel (1983) for monochromatic waves. We conjecture that this combination of beating between the bound and free superharmonic waves and modulation by their respective envelopes with each travelling at a different speed is the cause of the local peak in skewness and kurtosis near a depth transition reported in a series of previous papers reviewed in Trulsen et al. (2020).