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).

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 1983Massel , 1993Belibassakis & Athanassoulis 2002 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).
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π/(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).

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 h s , and the water depths can be deep (kh 1, with k the wavenumber), intermediate (kh = O(1)) or shallow (kh 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: where Φ(x, z, t) is the velocity potential and ζ(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: where g is the gravitational acceleration; a bottom boundary condition continuity of the potential and its horizontal derivative in the fluid exactly above the step (2.4a,b) and a no-flow boundary condition on the step wall

Stokes and multiple-scales expansions
In order to solve the boundary value problem (2.1)-(2.5), the unknown Φ and ζ are expressed as series solutions in the wave steepness = k 0 A (a so-called Stokes expansion), with k 0 and A denoting the characteristic wavenumber and wave amplitude, respectively: 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 , 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 ψ 0 = k 0 x 0 − ω 0 t 0 + μ 0 be the phase of the carrier wave, where ω 0 is the angular wave frequency, μ 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 = δx 0 and T = δt 0 are the slow scales and δ 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(δ) = O( )), 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 . 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) 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. Following Mei et al. (1989) and Massel (1983), we express the incoming wavepacket to leading order as

Description of the incoming wavepacket
which is valid for x 0, i.e. over the flat seabed to the left of the step. The superscript (mn, j) denotes the term of O( m δ j ) that is proportional to the harmonic exp(inψ 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(inψ 0 ) is understood). An analogous equation to (2.8) describes the free-surface elevation of the incoming wavepacket ζ I , and we proceed to express all the solutions in terms of the packet of its lowest-order term. Specifically, we assume 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) where κ m (0 < κ m 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.

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  (20,1) Match at x = 0 Match at  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.
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.
At first order in wave steepness, specifically O( δ 0 ), an incident wavepacket responds to an abrupt depth change by being reflected (Φ (11,0) R ) and transmitted (Φ (11,0) T ), complemented by the generation of evanescent waves (Φ (11,0) Ed on the deeper side and Φ (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( 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( 2 δ 0 )), which are proportional to exp(2iψ), and subharmonic bound waves (O( 2 δ 1 )), which are independent of the rapidly varying phase ψ 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.

First-order solutions (up to O( δ 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), Φ (1) is expressed as 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 Φ (11,1) does not contribute to the second-order solutions at O( 2 δ), but only to those at higher orders in bandwidth (see § 2.6). The details of the derivation of Φ (11,1) are nevertheless included in Appendix A for completeness. The linearised boundary value problem (2.1)-(2.5) yields as well as the evanescent wave coefficients R n and T m are complex, and I denotes the imaginary component. The coefficients R 0 , R n , T 0 and T m of the free waves at O( δ 0 ) are solved for numerically based on the boundary conditions at the step described by (2.4) and (2.5): 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 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 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: 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(−|ik n x|) or exp(−|ik m x|).
2.6. Second-order solutions (O( 2 )) The free-surface boundary conditions can be combined into one, which gives at O( 2 ) (e.g. Longuet-Higgins & Stewart 1964;McAllister et al. 2018) with a corresponding diagnostic equation for the surface elevation at O( 2 ) 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( 2 ) yields in which F 2 δ 0 and F 2 δ 1 can be further decomposed based on wave harmonics. A similar equation can be obtained for ζ (2) (not shown).
2.6.1. Superharmonic packets at O( 2 δ 0 ) Similar to Massel (1983), we seek solutions for the superharmonic packets at O( 2 δ 0 ) of the form 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 where the last two equations denote the dispersion relationships associated with the frequency 2ω 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 0s, 1, 2, . . .).

Subharmonic packets at O( 2 δ 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 2 δ 0 = 0. The solutions at O( δ) do not contribute to leading order, and we obtain the following for the subharmonic forcing at second order: where c 0s = ω 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: where κ 0 = −Ω/c g0 , κ 0s = −Ω/c g0s , Ω m (0 < δΩ m ω 0 ) is the maximum frequency of the packet resulting from the assumption of narrow bandwidth and (1 − tanh 2 k 0s h s ) + 1 . (2.29c) 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δ) 1) In the limit 1/(khδ) 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 u (20,1) where the non-dimensional coefficients B d and B s are given by When the carrier waves are additionally assumed to travel in deep water (i.e. k 0 h d 1 and k 0s h s 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/σ → 0, the order of the solutions in δ 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.
√ gh d on the deeper side and √ gh s on the shallower side. Assuming such free subharmonic packets can propagate in both directions, we seek solutions of the form ζ (20,1) (2.33b) The relationship between u and ζ is set by ∂ t u = −g∂ x ζ (cf. (4.1.3) in Mei et al. (1989)).
The coefficients B f R and B f T 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: where we note that u (1) (z = 0)η (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.35b) We note that the relations sign(B f T ) = −sign(B s ) and sign(B f R ) = −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 f T and B f T only depend on two non-dimensional parameters: k 0 h d and k 0 h s . We further explore these solutions and the underlying physics in the next section.

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: in which k 0 and ω 0 are the carrier wavenumber and angular frequency, respectively, is the group velocity on the deeper side, σ 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 = k 0 A 0 = 0.03 and the bandwidth parameter δ = 1/(k 0 σ 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. 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.

Generation of free packets: stage I versus stage III
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 ).
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 √ 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 √ gh d on the deeper side.

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 0 h 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  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 0 h d 2.0. The coefficient of reflection can reach a maximum of ∼ 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 ( 0.05π) and the reflected waves have a phase shift 0.2π when their amplitudes are ∼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π and −π. 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.

Behaviour near the abrupt depth transition: stage II
As noted in § 3.2, the transmitted superharmonic free wavepacket has a phase shift of approximately π 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 for any positive integer n with the first beat corresponding to n = 1, noting that arg(T 0 ) ≈ 0 and arg(T 2 ) ≈ −π. Taking 4σ x,s with σ x,s = σ 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 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) 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 λ 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 1.5. At least one beat can be expected for k 0s h s > 0.2 and δ s = 1/(k 0s σ x,s ) < 0.1, as shown in figure 6(b). As the group length increases (i.e. δ s decreases), more beats can be expected. When we eventually approach to the limit δ s → 0 (not shown here), denoting a uniform Stokes wave as studied and examined in Massel (1983), there are an infinite number of beats.

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).   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λ 0 and 80λ 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 1980Tayfun , 1986Mei 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.

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 915 A71-19  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).
where R 1i and T 1j ((i, j) = 0, 1, 2 . . .) are the unknown reflection and transmission coefficients of the free waves at this order (in analogous fashion to those at O( δ 0 )). Finally, at O( δ 1 ), the step boundary conditions (2.4) and (2.5) should be satisfied, which leads to simultaneous equations for the unknown coefficients R 1i and T 1j ((i, j) = 0, 1, 2 . . .). These can be solved for in a manner similar to the numerical method presented in Appendix B. We note that the solutions derived in this section can be readily checked by substituting the linear solutions up to O(δ) back into the linearised boundary value problem (2.1)-(2.5).
Appendix B. Numerical approach to obtain R n and T m As noted in § 2.5, the coefficients R n and T m are solved based on the step boundary conditions described by (2.13). We note that the following orthogonal properties of the hyperbolic (cosine) functions apply: in which h = h d (or h s ) and k i = k n (or k m ). Therefore, integrating the boundary conditions over the water column at x = 0 yields in which δ 0,0 = 1 and δ 0,n = 0 for n / = 0. Equation (B2) consists of N + M + 2 linear equations of N + M + 2 unknowns (i.e. R n and T m ), which can be solved numerically as a system of simultaneous equations.
Appendix C. Spatio-temporal dependence of the wavepackets Because the boundary conditions associated with the step (2.4) and (2.5) require evaluation of the solution at x = 0, information about the spatial dependence of the solution is lost at this stage. Letting the coefficients R 0 , T 0 , R n and T m capture magnitude and phase, the fact that the boundary conditions at the step (2.4) and (2.5) have to be satisfied for all time is captured by the time dependence of the packets at x = 0: (C1a-d) The spatial dependence of the packets can be obtained from the solvability condition (A2). Consequently, all (first-order) packets should vary as a function of X − c gi T, where c gi denotes the group velocity of the relevant packet, which is imaginary for evanescent waves. Inserting the linear potential at O( δ 0 ) into the boundary conditions at the step yields ik m T m cosh k j (z + h s ) cosh k j h s + ik 0s T 0 A T cosh k 0s (z + h s ) cosh k 0s h s for − h s < z < 0, (C2c) Applying orthogonality properties to (C2) (see Appendix B), the boundary conditions at the step can be rearranged in a matrix form as follows: where Y = [R 0 A R , R 1 A Ed,1 , R 2 A Ed,2 , . . . , R N A Ed,N , T 0 A T , T 1 A Es,1 , . . . , T M A Es,M ] T is the vector of unknowns, the matrix of coefficients C is composed of elements and K is a vector composed of elements given by The detailed expressions for k i h i (k j h j ) in (C4) and (C5) can be found in (B2) in Appendix B. Thus, we obtain for the unknowns The time dependence of Y and thus that of A R , A Ed,n , A T and A Es,m originates from the time dependence of A I . The coefficients are thus time independent and given by which is the same as (B2).