Long-Wave Penetration through a Laterally Periodic Continental Shelf

Abstract

The transformation of long waves—such as tsunamis and storm surges—evolving over a continental shelf is investigated. We approach this problem numerically using a pseudo-spectral method for a higher-order Euler formulation. Solitary waves and undular bores are considered as models for the long waves. The bathymetry possesses a periodic ridge-valley configuration in the alongshore direction which facilitates a means by which we may observe the effects of refraction, diffraction, focusing, and shoaling. In this scenario, the effects of wave focusing and shoaling enhance the wave amplitude and phase speed in the shallower regions of the domain. The combination of these effects leads to a wave pattern that is atypical of the usual behavior seen in linear shallow-water theory. A reciprocating behavior in the amplitude on the ridge and valley for the wave propagation causes wave radiation behind the leading waves, hence, the amplitude approaches a smaller asymptotic value than the equivalent case with no lateral variation. For an undular bore propagating in one dimension over a smooth step, we find that the water surface resolves into five different mean water levels. The physical mechanisms for this phenomenon are provided.

1. Introduction

There have been many studies to quantify the evolution of localized long waves during shoaling onto a continental shelf. For example, Synolakis [1] presented the transformation of solitary waves over plane beaches. Synolakis & Skjelbreia [2] proposed a “two-zone” shoaling model by performing laboratory experiments which show that at the commencement of the shoaling process (near the toe of the beach) that the rate of shoaling follows Green’s law which is applicable for small-amplitude shallow-water waves propagating over slowly-varying water depths. Green’s law can be written as $a\propto h^{-1/4}$ in which a is the wave amplitude and h is the local still water depth. After the solitary wave advances far enough towards the beach, the amplitude begins to coincide with a faster rate of shoaling ($a\propto h^{-1}$) corresponding to the adiabatic result by Miles [3]. Numerical results on solitary-wave shoaling generated by a fully-nonlinear wave model are presented by Grilli et al. [4] and Guyenne & Nicholls [5], which reveal the various rates of shoaling that may be expected for different beach slopes and wave amplitudes. A comprehensive numerical study of solitary-wave runup by Knowles & Yeh [6] reveals that the two-zone shoaling model is a special case among numerous shoaling behaviors depending on the length and slope of the beach, as well as the initial incident-wave amplitude. In some cases, the rate of wave shoaling can be slower than what is predicted by Green’s law.

The majority of the studies that have dealt with this problem approach it from the context of shoaling over a plane beach. However, some literature is available which includes the effects of geometric irregularities in the bathymetry. Chen et al. [7] numerically solved the extended Boussinesq equation [8] in the 2D horizontal plane for solitary waves in a nonuniform bathymetric domain. Solitary-wave-breaking over a circular shoal and run-up onto a conical island were examined and compared with results from large-scale laboratory experiments. Lynett et al. [9] performed laboratory experiments on solitary waves propagating over a shallow-water shelf. The continental shelf possessed a triangular shape which broadens in the nearshore direction. This bathymetry was modified by Lynett et al. [10] to contain a conical island located near the front of the triangular bathymetry. From these experiments, it was observed that the solitary wave would separate into four different bores. In the shallower regions, breaking would occur, generating jets around the conical island which led to an enhanced energy dissipation. Liu et al. [11] performed large-scale laboratory experiments on solitary waves propagating around a circular island. Results were compared with a numerical model of the shallow-water wave equations. The experiments provided evidence of wave trapping and the enhancing effects of slope on wave amplitude. Numerically, Yoshiki et al. [12] developed a depth-integrated, non-hydrostatic, shock-capturing model in spherical coordinates to reconstruct the 2009 Samoa tsunami. An analytic study was presented by Shimozono [13] in which solutions were developed for long-wave evolution in laterally converging bays. It was found in some cases that the wave behavior becomes two-dimensional due to wave refraction, and as a consequence the quasi-one-dimensional assumption breaks down.

Here we investigate the effects of alongshore periodic variations of continental shelves on long wave penetrations. This does not exactly represent a real-world condition, but we intend to focus on the physics that emerge from the presence of a 2D bathymetry.

Solitary waves are a convenient waveform used to study localized long waves because of their stable and permanent form in waters of uniform depth. Another relevant waveform that we shall consider are undular bores. An undular bore is a dispersive shockwave that manifests itself between two different asymptotic depths in shallow waters. The formation of an undular bore can emerge from an extremely long initial water-surface displacement such as a tsunami or storm surge. Furthermore, an undular bore is often generated by tides propagating up a river estuary. Similar to a fully turbulent bore, an undular bore is characterized by a sudden change in water depth which transitions from a supercritical to subcritical flow regime. It was found in laboratory experiments by Favre [14], that as long as the ratio between the upstream and downstream water depths is greater than unity but less than 1.28, an undular bore may form without any signs of energy dissipation. Peregrine [15] performed numerical experiments on undular bores in waters of uniform depth. He found that the initial growth rate of the undular bore is proportional to the square of its initial amplitude, and the calculations were compared to the experimental results found by Sandover & Taylor [16].

The purpose of our study is to explore the behaviors of a long wave (modeled by a solitary wave and an undular bore) while it climbs over a continental shelf with a lateral periodic variability. High-accuracy numerical experiments are carried out using the pseudo-spectral method. Quantitative analyses are performed which reveal the mechanisms responsible for the evolution of long waves over continental shelves.

2. Materials and Methods

Numerical Algorithm

We consider an irrotational flow field with a homogeneous, inviscid fluid in which the mathematical description is given by the Zakharov-Craig-Sulem formulation [17,18] of the full water-wave theory in the 3D domain:

$$\begin{array}{cc}{\mathsf{\Phi }}_{xx}+{\mathsf{\Phi }}_{yy}+{\mathsf{\Phi }}_{zz}=0\hfill & \mathrm{for} -1+\zeta \le z\le \eta ,\hfill \\ \nabla \mathsf{\Phi }·\nabla \zeta -{\mathsf{\Phi }}_z=0\hfill & \mathrm{at} z=-1+\zeta ,\hfill \\ \left.\begin{array}{c}{\displaystyle {\eta }_t+\nabla {\mathsf{\Phi }}^S{·\nabla \eta -(1+|\nabla \eta |}^2){\mathsf{\Phi }}_z=0}\hfill \\ {\displaystyle {\mathsf{\Phi }}_t^S+\eta +\frac{1}{2}|\nabla {\mathsf{\Phi }}^S{|}^2-\frac{1}{2}{(1+|\nabla \eta |}^2){\mathsf{\Phi }}_z^2=0}\hfill \end{array}\right\}\hfill & \mathrm{at} z=\eta ,\hfill \end{array}$$

(1)

where parameters have been scaled as follows:

$$\tilde{\mathsf{\Phi }}=c_0{h}_0\mathsf{\Phi }, \tilde{x}=h_{0}x, \tilde{y}=h_{0}y, \tilde{z}=h_{0}z, \tilde{\eta }=h_0\eta , \tilde{\zeta }=h_0\zeta , \tilde{t}=\frac{h_0}{c_0}t,$$

(2)

in which the tildes denote dimensional variables, $c_0=\sqrt{g{h}_0}$ is the linear shallow-water-wave phase speed, and $h_0$ is the still water depth. In this formulation, $\mathsf{\Phi }(\overrightarrow{x},z,t)$ is the velocity potential, $\zeta (\overrightarrow{x})$ is the deviation of the bottom-boundary from its mean elevation, $\eta (\overrightarrow{x},t)$ is the water-surface displacement from the equilibrium state, ${\mathsf{\Phi }}^S(\overrightarrow{x},t)\equiv \mathsf{\Phi }(\overrightarrow{x},z=\eta (\overrightarrow{x},t),t)$ is the velocity potential at the water surface, $\nabla =(\partial /\partial x,\partial /\partial y)$ represents the horizontal gradient, z points upward from the quiescent water level, and $\overrightarrow{x}=(x,y)$ is the horizontal position vector. From hereinafter, terms will be scaled consistently with (2). The velocity potential is expressed as an asymptotic series expansion:

$$\mathsf{\Phi }(\overrightarrow{x},z,t)=\sum _{m=1}^M{\varphi }^{(m)}(\overrightarrow{x},z,t),$$

(3)

where the superscript $(m)$ denotes that the term is on the order $O\sim ({\alpha }_0^m)$ in which ${\alpha }_0$ is a small parameter that represents the nonlinearity effect (${\alpha }_0=a_0/h_0$: $a_0$ is the wave-amplitude scale). Each perturbation term is expressed as a two-term expansion ${\varphi }^{(m)}=A^{(m)}+B^{(m)}$ in which:

$$A^{(m)}(\overrightarrow{x},z,t)=\sum _{p=-\infty }^{\infty }\sum _{q=-\infty }^{\infty }{A}_{p,q}^{(m)}(t)\frac{\cosh |{\overrightarrow{\kappa }}_{p,q}|(z+1)}{\cosh |{\overrightarrow{\kappa }}_{p,q}|}{e}^{i{\overrightarrow{\kappa }}_{p,q}·\overrightarrow{x}},$$

(4)

and

$$B^{(m)}(\overrightarrow{x},z,t)=\sum _{p=-\infty }^{\infty }\sum _{q=-\infty }^{\infty }{B}_{p,q}^{(m)}(t)\frac{\sinh |{\overrightarrow{\kappa }}_{p,q}|z}{|{\overrightarrow{\kappa }}_{p,q}|\cosh |{\overrightarrow{\kappa }}_{p,q}|}{e}^{i{\overrightarrow{\kappa }}_{p,q}·\overrightarrow{x}}.$$

(5)

The basis functions $A^{(m)}$ and $B^{(m)}$ satisfy Laplace’s equation in which (4) also satisfies the Neumann boundary condition at $z=-1$ and (5) satisfies the Dirichlet boundary conditon at $z=0$. For numerical considerations, the indices p and q of the basis functions are truncated to some large integers $N_x$ and $N_y$, respectively. In (4) and (5), ${\overrightarrow{\kappa }}_{p,q}=(\frac{p\pi }{L_x},\frac{q\pi }{L_y})$ is the wavenumber vector where $L_x$ and $L_y$ are the length and width of the domain, respectively. The magnitude of the wavenumber vector is $|{\overrightarrow{\kappa }}_{p,q}|=\sqrt{{(\frac{p\pi }{L_x})}^2+{(\frac{q\pi }{L_y})}^2}$, and $i=\sqrt{-1}$. A four-stage Runge-Kutta (RK4) method (e.g., [19]) is used to time-integrate the evolution equations numerically:

$$\begin{array}{cc}\hfill & {\eta }_t=-\nabla {\mathsf{\Phi }}^S{·\nabla \eta +(1+|\nabla \eta |}^2){\mathsf{\Phi }}_z,\hfill \\ \hfill & {\mathsf{\Phi }}_t^S=-\eta -\frac{1}{2}|\nabla {\mathsf{\Phi }}^S{|}^2+\frac{1}{2}{(1+|\nabla \eta |}^2){\mathsf{\Phi }}_z^2,\hfill \end{array}$$

(6)

on $z=\eta$. The horizontal gradient is calculated through the algebraic derivative by making use of the Fast Fourier Transform (FFT). In order to calculate the vertical velocity ${\mathsf{\Phi }}_z$ at the water surface, we must first calculate the degrees of freedom of the basis functions. In other words, we use the Fourier Transform to determine the coefficients $A_{p,q}^{(m)}(t)$ and $B_{p,q}^{(m)}(t)$. This is achieved by expanding the velocity potential at the mean water surface as well as at the bottom boundary. At the water surface:

$${\mathsf{\Phi }}^S(\overrightarrow{x},t)=\sum _{m=1}^M\sum _{k=0}^{M-m}\frac{{\eta }^k}{k!}\frac{{\partial }^k{\varphi }^{(m)}}{\partial z^k}{|}_{(\overrightarrow{x},z=0,t)},$$

(7)

and for example, if we set $M=5$, after expanding out the sum and equating like-order terms we have:

$$\begin{array}{c}{A}^{(1)}={\mathsf{\Phi }}^S,\hfill \\ A^{(2)}=-\eta {(A^{(1)}+B^{(1)})}_z,\hfill \\ A^{(3)}=-\eta {(A^{(2)}+B^{(2)})}_z-\frac{{\eta }^2}{2!}{A}_{zz}^{(1)},\hfill \\ A^{(4)}=-\eta {(A^{(3)}+B^{(3)})}_z-\frac{{\eta }^2}{2!}{A}_{zz}^{(2)}-\frac{{\eta }^3}{3!}{(A^{(1)}+B^{(1)})}_{zzz},\hfill \\ A^{(5)}=-\eta {(A^{(4)}+B^{(4)})}_z-\frac{{\eta }^2}{2!}{A}_{zz}^{(3)}-\frac{{\eta }^3}{3!}{(A^{(2)}+B^{(2)})}_{zzz}-\frac{{\eta }^4}{4!}{A}_{zzzz}^{(1)},\hfill \end{array}$$

(8)

all evaluated at $z=0$. The coefficients of (4) and (5) can then be found from the Fourier Transform once the $B^{(m)}$ terms are known. Therefore, we continue in a similar fashion by expanding the velocity potential about the bottom boundary:

$$\mathsf{\Phi }(\overrightarrow{x},z,t)=\sum _{m=1}^M\sum _{k=0}^{M-m}\frac{{(z+1)}^k}{k!}\frac{{\partial }^k{\varphi }^{(m)}}{\partial z^k}{|}_{(\overrightarrow{x},z=-1,t)}.$$

(9)

Substituting (9) into the kinematic boundary condition at the bed (second line of (1)) and equating like-order terms yields:

$$\begin{array}{c}{B}_z^{(1)}=0,\hfill \\ B_z^{(2)}=\nabla \zeta ·\nabla A^{(1)}-\zeta A_{zz}^{(1)},\hfill \\ B_z^{(3)}=\nabla \zeta ·\nabla (A^{(2)}+B^{(2)})-\zeta (A_{zz}^{(2)}+B_{zz}^{(2)}),\hfill \\ B_z^{(4)}=\nabla \zeta ·(\nabla A^{(3)}+\nabla B^{(3)}+\zeta \nabla B_z^{(2)}+\frac{{\zeta }^2}{2!}\nabla A_{zz}^{(1)})\hfill \\ -\zeta (A_{zz}^{(3)}+B_{zz}^{(3)}+\frac{\zeta }{2!}{B}_{zzz}^{(2)}+\frac{{\zeta }^2}{3!}{A}_{zzzz}^{(1)}),\hfill \\ B_z^{(5)}=\nabla \zeta ·(\nabla A^{(4)}+\nabla B^{(4)}+\zeta \nabla B_z^{(3)}+\frac{{\zeta }^2}{2!}\nabla A_{zz}^{(2)}+\frac{{\zeta }^2}{2!}\nabla B_{zz}^{(2)})\hfill \\ -\zeta (A_{zz}^{(4)}+B_{zz}^{(4)}+\frac{\zeta }{2!}{B}_{zzz}^{(3)}+\frac{{\zeta }^2}{3!}{A}_{zzzz}^{(2)}+\frac{{\zeta }^2}{3!}{B}_{zzzz}^{(2)}),\hfill \end{array}$$

(10)

all evaluated at $z=-1$. From (5), we see with the use of the Fourier Transform that the degrees of freedom $B_{p,q}^{(1)}(t)$ are all zero and therefore $B^{(1)}$ is zero for all t. Therefore $B^{(1)}$ may be omitted from (8) as well. From the recurrence relations established in (8) and (10) we can determine all of the degrees of freedom. After the degrees of freedom are determined, the vertical velocity at the water surface is now calculated:

$${\mathsf{\Phi }}_z(\overrightarrow{x},z=\eta ,t)=\sum _{m=1}^M\sum _{k=0}^{M-m}\frac{{\eta }^k}{k!}\frac{{\partial }^{k+1}{\varphi }^{(m)}}{\partial z^{k+1}}{|}_{(\overrightarrow{x},z=0,t)}.$$

(11)

The solution is completely determined once the initial conditions for $\eta$ and ${\mathsf{\Phi }}^S$ are prescribed. Note that the water-surface data must satisfy periodic boundary conditions to accommodate for the Fourier basis functions. This issue is resolved by mirroring the domain in the horizontal directions. The numerical method described here is based on the work by Dommermuth & Yue [20] and Gouin et al. [21]. The verification and validation of the present numerical model can be found in Knowles & Yeh [6].

3. Results

3.1. Solitary Wave Propagation over a Continental Shelf with Periodic Lateral Variation

Numerical results are presented for a solitary wave propagating over a continental shelf with a ridge-valley configuration as seen in the planview image presented in Figure 1. The bathymetry is given by:

$$\zeta (x,y)=\frac{{\delta }_0}{2}(1+\tanh {\sigma }_x[x+b_0\cos {\sigma }_{y}y-l_0]),$$

(12)

where the height of the step is ${\delta }_0=0.5$, the steepness in the $x-$direction is ${\sigma }_x=0.01$, and the amplitude of the lateral perturbation in the bathymetry is $b_0=50$. In order to ensure periodicity in the $y-$direction, we set ${\sigma }_y=\frac{2\pi }{L_y}=6.14\times {10}^{-2}$: namely the breadth between valleys is $L_y=102.2$. The bathymetry is translated to the right of the origin by an amount of $l_0=480$ so that the initial wave form is not disturbed by the presence of the continental shelf. Note that our coordinate system is oriented such that x is in the dominant direction of wave propagation where the origin $(x=0)$ gives the initial location of the crest of the solitary wave, and y is the spatial coordinate in the lateral direction. The bathymetry is symmetric about the origin $(y=0)$ on account of the spatial periodicity in the $y-$direction. Because of this periodicity, from hereinafter, the numerical results are presented only in the span $|y|\le 51.1$, and we call the bathymetry presented in (12) the 2D case, while the bathymetry with no lateral variation will be referred to as the 1D case.

The spatial resolutions in the $x-$ and $y-$directions are chosen to be $\Delta x=0.37$ and $\Delta y=0.1$, respectively. A time step of $\Delta t=0.16$ is used in this simulation and we truncate the asymptotic series (3) to $M=5$ terms. We include $N_x=2^{12}$ and $N_y=2^{10}$ modes in the basis functions (4) and (5).

The initial conditions are prescribed by the third-order solitary-wave solution by Grimshaw [22] where the initial water-surface displacement $\eta$ and velocity potential ${\mathsf{\Phi }}^S$ are given by:

$$\begin{array}{cc}\hfill \eta & ={\alpha }_0{S}^2-\frac{3}{4}{\alpha }_0^2(S^2-S^4)+{\alpha }_0^3\left(\frac{5}{8}{S}^2-\frac{151}{80}{S}^4+\frac{101}{80}{S}^6\right),\hfill \\ \hfill {\mathsf{\Phi }}^S& =2\sqrt{\frac{{\alpha }_0}{3}}\{T+{\alpha }_0\left[\frac{5}{24}T-\frac{1}{3}{S}^{2}T+\frac{3}{4}{(1+\eta )}^2{S}^{2}T\right]+\hfill \\ & {\alpha }_0^2[-\frac{1257}{3200}T+\frac{9}{200}{S}^{2}T+\frac{6}{25}{S}^{4}T+{(1+\eta )}^2\left(-\frac{9}{32}{S}^{2}T-\frac{3}{2}{S}^{4}T\right)+\hfill \\ \hfill & {(1+\eta )}^4(-\frac{3}{16}{S}^{2}T+\frac{9}{16}{S}^{4}T)\left]\right\}.\hfill \end{array}$$

(13)

where S and T are:

$$S=\mathrm{sech}{\kappa }_0(x-x_0), T=\tanh {\kappa }_0(x-x_0),$$

(14)

in which ${\kappa }_0=\frac{\sqrt{3{\alpha }_0}}{2}(1-\frac{5}{8}{\alpha }_0+\frac{71}{128}{\alpha }_0^2)$. In this experiment, the initial wave amplitude is ${\alpha }_0=0.1$. Recall that ${\alpha }_0$ is the wave amplitude normalized by the offshore still water depth and is therefore a measure of the initial nonlinearity of the wave.

In Figure 2, we present snapshots of the water-surface displacement for the solitary wave evolving over the bathymetry specified by (12). In Figure 2a, we present the initial offshore solitary-wave conditions. Then at time $t=470$ (Figure 2b), the effects of refraction begin to emerge from the presence of the bathymetry, creating a phase difference between the waveform at the ridge ($y=0$) and at the valley ($y=\pm 51.1$). Note that soon after refraction occurs, a stem-wave formation appears in the vicinity of $y=0$, which is similar in form to the wave pattern created by obliquely interacting solitary waves in a uniform depth: see for example, Yeh & Li [23] and Kodama & Yeh [24]. Due to the formation of “radiating” waves from the stem, energy is allowed to leak out laterally which is then displaced behind the lead waveform.

At time $t=658$ (Figure 2f), we observe a negative difference in spatial phase lag ($\Delta <0$) even though the local water depth is shallower at $y=0$ than at $y=\pm 51.1$. Note that the spatial phase difference is defined as $\Delta =x_v-x_c$, where $x_v$ gives the $x-$coordinate of the lead crest of the solitary wave at $y=\pm 51.1$, and $x_c$ gives the $x-$coordinate of the lead crest at $y=0$. This negative phase lag is also observed in Figure 2g,h. While the negative phase lag is present, energy is allowed to travel outwards in the lateral directions, transverse to the direction of wave propagation, which then accumulates at $y=\pm 51.1$, at the valleys of the bathymetry. This then commences in an increase in amplitude and thus phase speed at the valleys, which leads to a positive spatial phase difference ($\Delta >0$) as can be observed at $t=940$ (Figure 2i). In Figure 2j,k we continue to see a reciprocating behavior in the spatial phase lag where the middle of the waveform bends forwards and backwards in time. However, as the wave climbs up the shelf, the waveform is tending towards an asymptotic state where the solitary wave is recovered and the crest of the wave becomes straight and parallel to the shoreline. Note that the waveform at time $t=1409$ (Figure 2k) is of greater amplitude, but of narrower width, than the initial solitary wave. The total mechanical energy of the initial solitary wave is compared with the lead solitary wave at the end of the numerical simulation, and the ratio of the final-stage energy to the initial energy was found to be 0.74. This is the amount of energy reduction in the leading wave due to dispersion/radiation effects which took place over the wavy bathymetry.

In Figure 3a, we present a portion of the domain along the ridge ($y=0$) at time $t=2036$. The waveform over the plateau region possesses a lead solitary wave and trailing waves of smaller amplitude. In order to assess if these are truly solitary waves, we compare the computed numerical water-surface displacement with the higher-order solution (13) by Grimshaw [22]. An overlay presented in Figure 3 shows that the lead and intermediate waves are solitary waves. However, the third wave, presented in Figure 3d, does not appear to be a clean solitary wave.

In Figure 4, we present the amplitude variation along the ridge $(y=0)$ and the valley $(y=\pm 51.1)$. In addition, the result from the 2D case is compared with the 1D wave-propagation case which will be discussed later in Section 3.1.1. The amplitude at $y=0$ starts to increase at the toe of the slope ($x\le 250$) and reaches a maximum of $a=0.274$ at $x=556$. The amplitude at $y=\pm 51.1$ then attains a minimum value of $a=0.038$ at $x=651$. We observe a “saw-tooth” pattern where the amplitude at $y=0$ oscillates out of phase with the amplitude at $y=\pm 51.1$. The amplitude appears to be approaching an asymptotic mean state of $a\approx 0.156$ for the 2D case, whereas the amplitude for the 1D case is approaching a larger value of $a\approx 0.183$. This smaller amplitude for the 2D case is due to energy dispersion effects resulting from the creation of radiating waves in the vicinity of $y=0$.

The spatial phase difference $\Delta =x_v-x_c$ is presented in Figure 5. Again, we see an oscillatory behavior in the phase lag which also appears to be tending towards an asymptotic state: $\Delta \to 0$ as $x\to \infty$. The maximum positive value of phase lag is $\Delta =6.18$ at the location $x=495$. The most negative value of phase lag is $\Delta =-7.59$ at the location $x=696$. Note that, the greatest absolute difference in phase lag occurs when the waveform in the shallower region of the domain is moving faster than the waveform in the deeper portion of the domain.

3.1.1. Solitary Wave Propagation over a Continental Shelf with no Lateral Variation

For comparison, we run a numerical simulation for the 1D wave-propagation case. The initial wave conditions are given by equation (13) where the amplitude is ${\alpha }_0=0.1$ and the initial location of the crest is $x_0=0$. The bathymetry with no lateral variation is now given by:

$$\zeta (x)=\frac{{\delta }_0}{2}(1+\tanh {\sigma }_x[x+b_0-l_0]).$$

(15)

The parameters for the bathymetry (15) are ${\delta }_0=0.5$, ${\sigma }_x=0.01$, $b_0=50$, and $l_0=480$. Note that this is the same bathymetry as for the 2D wave-propgation case (12) if we had chosen $y=0$. A temporal and spatial resolution of $\Delta t=0.016$ and $\Delta x=0.07$ are used, respectively. We include $M=5$ terms in the perturbation expansion and truncate the basis functions to $N_x=2^{15}$ modes.

In Figure 6, we provide snapshots of the free-surface displacement along with the corresponding bathymetry at various times. Figure 6a shows the initial condition of the waveform at the offshore location. At time $t=470$ (Figure 6b), the solitary wave has encountered the increasing bed elevation and the water-surface displacement loses its symmetry. The amplitude increases and the solitary wave becomes narrower. In Figure 6c, we observe the presence of a trailing wave resulting from fissioning that occured after the solitary wave climbed onto the plateau region of the bathymetry. The solitary waveform is then recovered as demonstrated in Figure 6d with a larger amplitude and narrower width than the initial solitary wave.

The waveform is presented in Figure 7 at time $t=2036$. We see a similar pattern in the water-surface displacement along $y=0$ for the 2D case seen in Figure 3. However, for the 2D case, there are extra trailing waves. As we stated earlier, the amplitude of the lead solitary wave is larger for the 1D case than for the 2D case. Furthermore, in Figure 7, we overlay the numerical result with the higher-order solitary-wave solution (13) by Grimshaw [22] in order to check if the trailing and leading waves are solitary waves. It appears that the first two leading waves, as seen in Figure 7b,c are solitary waves, but the trailing wave, as seen in Figure 7d is similar to a solitary wave, but not quite the exact same form. Therefore, just as in the 2D case, it appears that there are two solitary waves that emerge from the fissioning process. The numerical result for the number of solitons which emerge is compared with the analytical prediction presented in (Johnson [25] Section 3.4.4):

$$1-{\delta }_0={\{\frac{1}{2}(N)(N+1)\}}^{-4/9},$$

(16)

which gives a value of $N=2.62$ for ${\delta }_0=0.5$. Our numerical results support the analytical prediction (16). Note that, in the present case we examined, the analytical prediction coincides with the numerical result even for the 2D bathymetry.

3.2. Undular Bore Propagation over a Continental Shelf with Periodic Lateral Variation

We now examine an undular bore propagating over a bathymetry representing a continental shelf with a laterally-periodic ridge-valley configuration, similar to what is seen in Figure 1. The initial wave condition is prescribed by a smooth step-shaped displacement with zero velocity. The initial water-surface displacement and velocity potential are given, respectively, by:

$$\eta =\frac{{\alpha }_0}{2}(1-\tanh {\kappa }_{0}x), {\mathsf{\Phi }}^S=0,$$

(17)

where ${\alpha }_0=0.1$ and ${\kappa }_0=1.0$. The bathymetry is prescribed by (12) where the parameters are given by ${\delta }_0=0.5$, ${\sigma }_x=0.01$, ${\sigma }_y=6.14\times {10}^{-2}$, and $b_0=50$. In this simulation we choose $l_0=980$ in order to give the initial wave condition ample time to evolve into an undular bore before encountering the continental shelf. We set a time step at $\Delta t=0.16$ and the spatial resolutions in the $x-$ and $y-$directions are $\Delta x=0.37$ and $\Delta y=0.10$, respectively. $M=5$ terms are included in the perturpation series, and the modes in the $x-$ and $y-$directions are truncated to $N_x=2^{13}$ and $N_y=2^{10}$, respectively.

Figure 8 shows snapshots of the plan view of the water-surface displacement along with the corresponding bathymetry. In Figure 8a, we see the initial offshore water-surface displacement at $t=0$. The smooth step then evolves into an undular bore in the offshore region as shown in Figure 8b. At time $t=940$ (Figure 8c), the undular bore encounters the toe of the continental slope and refracts about the crest of the ridge. In Figure 8d,e, the lead crest of the undular bore begins to straighten out while “radiating” waves emerge from the ends of the stem; this is similar to the pattern observed in the solitary wave case, as discussed in Section 3.1. At time $t=1096$ (Figure 8f), it appears that the center of the waveform begins to bow forward which continues in Figure 8g,h. At time $t=1253$, (Figure 8i), the amplitude at the valley ($y=\pm 51.1$) becomes greater than at the ridge. In Figure 8j, solitons disintegrate from the front of the undular bore, resulting in a total of five solitary waves that have emerged at time $t=2193$ as shown in Figure 8k.

Figure 9 shows the spatial evolution of the lead-wave amplitude of the 2D undular bore against location. The amplitude at the ridge ($y=0$) is represented by a solid line and the amplitude at the valley ($y=\pm 51.1$) is represented by a dashed line. The dotted line represents the amplitude of the lead crest of the undular bore for the 1D case which is discussed later in Section 3.2.1. The amplitude at $x=0$ is $a_s=0.049$, which represents the height immediately after the intial water-surface displacement is set in motion. This initial height is in agreement with the prediction given by the method of characteristics for the nonlinear-nondispersive shallow-water-wave equation for a dam-break initial condition:

$$a_s=\frac{{(1+\sqrt{{\alpha }_0+1})}^2}{4}-1.$$

(18)

Soon after, the amplitude of the lead wave of the undular bore slowly increases until reaching $x\approx 875$. Once the undular bore encounters the sloping bed, we see that the amplitude at the ridge ($y=0$) quickly increases, reaching a maximum value of $a=0.266$ at $x=1044$, and the amplitude in the valley $(y=\pm 51.1)$ decreases, reaching a minimum value of $a=0.035$ at $x=1160$. The amplitude then behaves in a reciprocating fashion, creating a saw-tooth pattern. The amplitude at $y=0$ is roughly out of phase with the amplitude at $y=\pm 51.1$. It appears though that the amplitude may approach an asymptotic state given a sufficient amount of time. For the 1D case that is represented by the dotted line in Figure 9, the amplitude increases monotonically until reaching the asymptotic value of $a=0.17$ after reaching the plateau. It appears that the amplitudes for the 1D and 2D cases become similar in the limit as $x\to \infty$. Recall that the amplitudes for the 1D and 2D cases of the solitary waves are different ($a\approx 0.183$ for the 1D case and $a\approx 0.156$ for the 2D case) as $x\to \infty$.

The spatial phase lag $\Delta =x_v-x_c$ is presented in Figure 10 for the 2D undular bore case. The phase lag remains at $\Delta =0$ until encountering the toe of the bed slope at $x\approx 875$. At this point, a positive phase lag is created from the undular bore refracting about the crest of the continental shelf. A maximum value of $\Delta =6.62$ is attained at $x=1012$. However, the spatial phase lag then sharply decreases to a minimum value of $\Delta =-6.46$ at $x=1197$. It appears that the phase lag is approaching an asymptotic value such that $\Delta \to 0$ as $x\to \infty$. Note that, unlike the solitary wave case (Figure 5), the magnitudes of the maximum and minimum phase lags are similar for the undular bore.

3.2.1. Undular Bore Propagation over a Continental Shelf with no Lateral Variation

For comparison, we perform numerical experiments where an undular bore propagates over a smooth step into shallower waters without lateral variation (see the sketch in Figure 11). The same initial conditions along $y=0$ as (17) is applied. In this numerical experiment, we set $\Delta t=0.016$, $\Delta x=0.13$, $M=5$, and $N_x=2^{16}$.

Snapshots of the entire domain are presented in Figure 12. The initial conditions are shown in Figure 12a. The smooth step then evolves into an undular bore and a dispersive rarefaction wave with an expanding mean water depth, which separates the two ends of the waveform. The mean water depth is equivalent to the height predicted by Equation (18). Once encountering the toe of the continental slope, the undular bore shoals, and a careful observation in Figure 12c reveals the occurrence of energy reflection at time $t=1096$. Note that, at the beginning of the simulation, there are only two mean water depths. However, after the undular bore climbs up and passes over the sloping region, five distinct mean water depths can be identified at $t=3132$ as shown in Figure 13. We may roughly partition these depths into the following regions: (1): $\overline{\eta }=0.1$ for $x<-3300$; (2): $\overline{\eta }=0.049$ for $-3000<x<-1600$; (3): $\overline{\eta }=0.058$ for $-800<x<750$; (4): $\overline{\eta }=0.055$ for $1200<x<2200$; and (5): $\overline{\eta }=0$ for $x>2800$. Regions (1), (2), and (5) are relatively straightforward to interpret. However, regions (3) and (4) require some explanations. Within region (4), the velocity field reaches an essentially steady flow condition. This is supported by Figure 14 which displays the depth-averaged horizontal velocity $\overline{u}$ plotted against time at locations $x=400$ and $x=1400$. After about $t=2500$, the velocity appears to maintain the steady state. Note that the depth-averaged horizontal velocity is defined as

$$\overline{u}=\frac{1}{H}{\int }_{-1+\zeta }^{\eta }u\mathrm{d}z,$$

(19)

in which u is the horizontal fluid velocity and $H=\eta +1-\zeta$ is the local total water depth. The fluid velocities are normalized by $c_0=\sqrt{g{h}_0}$. The vertical component of the velocity is essentially nil. For an inviscid and homogeneous fluid, the total energy head (total mechanical energy plus pressure work per unit fluid weight) must be conserved. Therefore:

$$E_3=E_4+\Delta z,$$

(20)

where $E_3$ and $E_4$ are the specific energy before and after the bed slope, respectively, which are defined as:

$$E_3=\frac{{\overline{u}}_3^2}{2}+H_3,E_4=\frac{{\overline{u}}_4^2}{2}+H_4,$$

(21)

where the subscripts 3 and 4 refer to the values in the regions 3 and 4, respectively, (see Figure 13), and H is the total water depth. Following this notation, the height of the step can be given by the difference in the bed elevations: $\Delta z={\zeta }_4-{\zeta }_3$. From Equations (20) and (21), the theoretical total water depth at region 4 can be computed and found to be 0.556. This is calculated from the numerical data attained at $t=3132$ with ${\overline{u}}_3=4.042\times {10}^{-2}$, ${\overline{u}}_4=7.67\times {10}^{-2}$, $H_3=1.058$, and $\Delta z=0.50$. The actual value attained from the numerical experiment was $H_4=0.556$, resulting in no error up to the first three decimal places.

The analysis between regions (3) and (4) indicates that the flow is subcritical relative to the laboratory coordinate, and a standard analysis of the specific energy in the area of open-channel flows (see e.g., Henderson [26]) can be applicable here for the quasi-steady flow condition.

4. Discussion

Numerical experiments are performed for a solitary wave and an undular bore evolving over a continental shelf containing a lateral variability. The repetitious ridge-and-valley pattern of the continental shelf causes the incident wave to deform due to the effects of refraction, diffraction, and shoaling. The local wave focusing produces the formation of a stem-like wave pattern that results from the spatially merging wave. The stem-like formation at the merging location emanates “radiating” trailing oblique waves from the edges of the stem. These local interactions on the ridge are similar to the interaction of V-shaped solitary waves in a uniform depth [27], which also produce trailing wave radiation. These radiating waves are the source of additional wave-energy dispersion. Consequently, the energy of the solitary wave climbing over the 2D continental slope is less than the comparable 1D case (see Figure 4). On the other hand, the energy loss in the leading wave of the undular bore is substantially smaller than the case of the solitary wave (Figure 9). This is because the wave emanating from the stem interacts with the subsequent wave of the undular bore; the radiated wave cannot freely disperse behind the leading wave. The accumulation in energy of the obliquely radiating wave tends to maintain the leading-wave amplitude of the undular bore. The foregoing explanation can be supported by the comparison of our numerical results shown in Figure 2f for the solitary wave and Figure 8f for the undular bore.

After the solitary wave reaches the plateau, soliton fissioning occurs, which results in the formation of two distinct solitary waves in the simulation. The numerical results of the soliton fission for solitary waves propagating over both the 1D and 2D cases agree with the analytical prediction given by Equation (16). After a sufficient amount of time, the solitary wave recovers its characteristic form, approaching a steadily propagating state. It is found that the amplitude of the 1D propagation case approaches a greater asymptotic value than for the 2D propagation case. This difference in amplitude appears to be a result of extra energy dispersion caused by the formation of the obliquely radiating waves, which is allowed to occur in the 2D case.

Results were presented on an undular bore evolving from dam-break initial conditions and then propagating over a smooth depth transition from deep to shallow water. Note that this initial condition resembles an ideal model of a co-seismic tsunami source. After the wave reaches the plateau region, multiple (in this case five) discrete solitary waves appear to disintegrate from the leading portion of the undular bore. This particular result has potential applications in engineering. For instance, after a tsunami propagates over a continental slope, our results demonstrate the possibility for multiple localized waves of significant amplitude to advance over a broad continental shelf and ultimately inundate the shore with a series of waves.

The present results demonstrate the uncertainty in predicting the critical locations of the greatest tsunami amplitude. From linear theory, we anticipate that the greatest amplitude should be observed along the ridge ($y=0$). However, our numerical experiments reveal that there is a laterally reciprocating variability in the amplification which occurs as a consequence of the nonlinear dependence of the phase speed on the amplitude.

The mean water surface of the entire domain was also examined for an undular bore propagating over the 1D continental shelf. Initially, we see that the water surface only contains two mean depths, but after the undular bore propagates well beyond the continental shelf, we find a total of five mean water depths over different segments of the domain. It appears that specific energy (a standard parameter used to analyze 1D steady open-channel flows) can be used to accurately predict the height of these mean water depths before and after the continental shelf.

5. Conclusions

Many previous studies on tsunami propagation over a continental slope assessed the problem for a 1D shelf configuration: i.e., no lateral variation. Here we examined long-wave (tsunami) transformations over a continental slope with a periodic lateral variation, considering two incident wave models: a solitary wave and an undular bore. The lateral variation causes the incident wave to refract, diffract and shoal along the trajectory of wave propagation. When the wave initially encounters the continental slope, shoaling occurs, placing the location of maximum amplitude over the ridge. Soon after, the location of maximum amplitude shifts to the valley due to the effects of wave nonlinearity and phase velocity in the 2D field. This transfer of energy then periodically shifts back and forth between the ridge and the valley. This process induces extra wave energy dispersion, and affects the energy distribution on the continental shelf after the incident wave climbs up the slope. Our findings demonstrate the importance in consideration of lateral variability of the continental slope for the evaluation of long-wave behaviors and consequences when it climbs up the continental slope.

The behaviors of localized long waves interacting and shoaling over a 2D bathymetry are primarily discussed in the context of length scales on the order of a continental shelf, where the corresponding waveform is a tsunami. None the less, the findings should be applicable to smaller scale nearshore waves interacting with the bathymetry as well.