2. The role of the sediment layer

The role of a porous sea bed in attenuating hydro-acoustic waves has been investigated by Chierici et al. (Reference Chierici, Pignagnoli and Embriaco2010). They proposed two theoretical solutions for constant water depth, the first based on compressible sea water and an incompressible sediment layer using Darcy equations. The second model is based loosely on the work of Buckingham (Reference Buckingham1997), who showed that the effect of intergranular stresses in an unconsolidated sediment plays the role of an apparent viscous dissipation added to the dissipation associated with the pore water percolation; this model serves as the basis for the derivation in § 3. Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013) have investigated the role of an elastic sea bed on progressive waves and found that the first acoustic mode ( $n=1$ ) is the dominant component of the hydro-acoustic wavefield. Here, we consider the interaction of a train of hydro-acoustic waves in a water column of depth $h(x,y,t)$ overlying a sediment layer of thickness $a(x,y,t)$ , with $h_{s}=h+a$ , where $h_{s}(x,y,t)$ is the total depth. The vertical coordinate, $z$ , is measured positively upwards from the undisturbed free surface at $z=0$ , and $x$ and $y$ denote horizontal Cartesian coordinates as shown in figure 1. The sediment layer causes the damping of hydro-acoustic waves, lowering the whole energy spectrum and shifting the expected frequency peaks towards lower values (Chierici et al. Reference Chierici, Pignagnoli and Embriaco2010). In the absence of viscous behaviour of the sea bottom, the dominant frequency range in the wave spectrum can be expressed by a discrete set of normal frequencies $f_{n}$ given by

(2.1) $$\begin{eqnarray}f_{n}=(2n-1)\frac{c}{4h},\quad n=1,2,3,\dots\end{eqnarray}$$

where $c$ is the sound speed in water (approximately $1500~\text{m}~\text{s}^{-1}$ ). Introducing the underlying sediment layer acting together with the water column causes lowering of the spectrum peaks determined from the following transcendental equation (Nosov et al. Reference Nosov, Kolesov, Denisova, Alekseev and Levin2007):

(2.2) $$\begin{eqnarray}\tan \left[\frac{2{\rm\pi}{\it\gamma}_{n}h}{c}\right]\tan \left[\frac{2{\rm\pi}{\it\gamma}_{n}a}{c_{s}}\right]=\frac{{\it\rho}_{s}c_{s}}{{\it\rho}c},\end{eqnarray}$$

where $c_{s}$ is the sound speed in sediment, ${\it\rho}$ and ${\it\rho}_{s}$ are the densities of water and sediment respectively, and ${\it\gamma}_{n}$ denotes the normal mode frequencies for the damped case. Note that in the case of $a=0$ , the set of normal modes described by (2.1) and (2.2) coincide.

Figure 1. Schematic view of fluid domain.

Figure 2. Bottom-pressure records at a distance of $x=96~\text{km}$ from the epicentre. Results of a sample computation carried out using a 3D flow solver in a constant depth, $h=2200~\text{m}$ , $a=1000~\text{m}$ , $c=1500~\text{m}~\text{s}^{-1}$ , $c_{s}=2000~\text{m}~\text{s}^{-1}$ , ${\it\rho}=1000~\text{kg}~\text{m}^{-3}$ , ${\it\rho}_{s}=1850~\text{kg}~\text{m}^{-3}$ for a unit source area with semi-length $b=112~\text{km}$ and rise time ${\it\tau}=1~\text{s}$ , showing: (a) time series of a one-layer compressible water model (light grey) and a coupled model with compressible water and inviscid compressible sediment (black), ${\it\mu}=0$ ; (b) the corresponding frequency spectrum to (a); (c) time series of a one-layer compressible water model with the partially reflecting boundary condition (2.4) at the bottom (light grey) and a coupled model with compressible water and viscous compressible sediment, ${\it\mu}_{s}=2\times 10^{8}~\text{Pa}~\text{s}$ (black); (d) the corresponding frequency spectrum to (c). The vertical dashed lines in (b,d) indicate the frequency peaks calculated by (2.1) and (2.2) in light grey and black respectively.

A sample computation is carried out using a full 3D solver in constant depth. The governing equations within layers and boundary conditions at free surface, interface and bottom for sample computation are described in § 3. We use $h=2200~\text{m}$ , $a=1000~\text{m}$ , $c=1500~\text{m}~\text{s}^{-1}$ , $c_{s}=2000~\text{m}~\text{s}^{-1}$ , ${\it\rho}=1000~\text{kg}~\text{m}^{-3}$ and ${\it\rho}_{s}=1850~\text{kg}~\text{m}^{-3}$ . The other parameters are for a unit sudden elevation of source area with semi-length $b=112~\text{km}$ and rise time ${\it\tau}=1~\text{s}$ . The transient sea bed velocity, $h_{t}$ , with a residual displacement $h_{0}$ , is a trigonometric function expressed by

(2.3) $$\begin{eqnarray}h_{t}=\frac{h_{0}}{2}\left[1-\cos \left(\frac{2{\rm\pi}(t-t_{0})}{{\it\tau}}\right)\right][H(t-t_{0})-H(t-t_{0}-{\it\tau})],\end{eqnarray}$$

where $H(t)$ is the Heaviside step function. The results are depicted in figure 2, which show the bottom pressure $P$ and the corresponding frequency spectrum $\tilde{P}$ at 96 km from the epicentre. In figure 2(a), the light grey line shows the time series of bottom pressure from a one-layer compressible water model while the black line shows the results from a model coupling two layers of compressible water and sediment without a damping term ( ${\it\mu}_{s}=0$ ) in the lower layer. In figure 2(b), for the one-layer model, the frequency peaks coincide with the cutoff frequencies for an ideal impermeable bottom ( $f_{1}=0.17$ and $f_{2}=0.51~\text{Hz},\dots$ ) identified by (2.1). For the case of the coupled model, the spectrum is peaked at ${\it\gamma}_{1}=0.15$ and ${\it\gamma}_{1}=0.407~\text{Hz}$ , representing cutoff frequencies evaluated by (2.2). In order to consider the damping behaviour of the underlying sediment layer, the bulk viscosity (ranging from $10^{6}$ up to $10^{20}~\text{Pa}~\text{s}$ (Van Keken et al. Reference Van Keken, Spiers, Van den Berg and Muyzert1993; Kimura Reference Kimura2006)) is fixed at $2\times 10^{8}~\text{Pa}~\text{s}$ . The model results are compared with a one-layer model with an additional partially reflecting boundary condition at the bottom defined by (2.4)

(2.4) $$\begin{eqnarray}{\it\Phi}_{z}=\frac{1}{c}\frac{K_{r}-1}{K_{r}+1}{\it\Phi}_{t},\end{eqnarray}$$

where $K_{r}$ is the amplitude reflection coefficient (Brekhovskikh, Lysanov & Lysanov Reference Brekhovskikh, Lysanov and Lysanov2003),

(2.5) $$\begin{eqnarray}K_{r}=\frac{{\it\rho}_{s}c_{p}-{\it\rho}c}{{\it\rho}_{s}c_{p}+{\it\rho}c},\end{eqnarray}$$

and $c_{p}$ is the P-wave speed in the bottom. This boundary condition causes decay of hydro-acoustic waves peaking at the cutoff frequencies ( $f_{n}$ ) defined by (2.1). In figure 2(c), the light grey line shows the time series of bottom pressure from the one-layer compressible model with a damping coefficient. Here, $c_{p}=8000~\text{m}~\text{s}^{-1}$ is selected for computations. Thus, the reflection coefficient takes the value of $K_{r}\approx 0.816$ . It can be seen from the plot that the oscillation lasts for 80 s after the earthquake. The lower time series is for a two-layer model with a dissipation term in the sediment (black). The attenuation of hydro-acoustic waves lasts for 500 s. The corresponding frequency spectra are shown in 2(d). Lowering of the whole energy spectrum in comparison with panel 2(b) is distinguishable. The comparison shown in figure 2 justifies the mismatch between the calculated spectral peaks for an impermeable bottom, $f_{n}$ , and the dominant frequencies, ${\it\gamma}_{n}$ , observed during the Tokachi Oki 2003 event (Nosov et al. Reference Nosov, Kolesov, Denisova, Alekseev and Levin2007).

3. The mild-slope equation for damped hydro-acoustic waves

We develop a mild-slope equation based on the eigenfunction structure for the problem with constant layer depths $h$ and $a=h_{s}-h$ and with no lower layer damping. We review the governing equations and boundary conditions here and then derive the damped, two-layer mild-slope model using an approach described in Silva, Salles & Govaere (Reference Silva, Salles and Govaere2003).

3.1. Governing equations

The linearized weakly compressible wave equation governing the fluid potential ${\it\Phi}(x,y,z,t)$ in the water layer is given by

(3.1a,b ) $$\begin{eqnarray}L({\it\Phi})={\it\Phi}_{tt}-c^{2}{\rm\nabla}^{2}{\it\Phi}=0;\quad -h+{\it\eta}_{2}(x,y,t)\leqslant z\leqslant {\it\eta}_{1}(x,y,t),\end{eqnarray}$$

where ${\rm\nabla}^{2}$ is the Laplacian in 3D and subscripts on dependent variables denote partial derivatives. The interfacial displacements ${\it\eta}_{1}$ and ${\it\eta}_{2}$ represent the response of the free surface and the layer interface to hydro-acoustic disturbances. Following Buckingham (Reference Buckingham1997) and Chierici et al. (Reference Chierici, Pignagnoli and Embriaco2010), the weakly compressible wave equation for the fluid potential $\mathscr{Q}(x,y,z,t)$ in the viscous sediment layer is given by

(3.2a,b ) $$\begin{eqnarray}L_{s}(\mathscr{Q})=\mathscr{Q}_{tt}-c_{s}^{2}{\rm\nabla}^{2}\mathscr{Q}-2{\it\nu}_{s}({\rm\nabla}^{2}\mathscr{Q})_{t}=0;\quad -h_{s}\leqslant z\leqslant -h+{\it\eta}_{2}(x,y,t)\end{eqnarray}$$

with apparent sediment kinematic viscosity ${\it\nu}_{s}$ . The boundary conditions at the free surface and bottom are given by

(3.3) $$\begin{eqnarray}\left\{\begin{array}{@{}c@{}}{\it\Phi}_{tt}+g{\it\Phi}_{z}=0\quad \text{at}~z=0\\ \mathscr{Q}_{z}+\boldsymbol{{\rm\nabla}}_{h}h_{s}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{h}\mathscr{Q}+h_{s,t}=0\quad \text{at}~z=-h_{s}\end{array}\right.\end{eqnarray}$$

where $\boldsymbol{{\rm\nabla}}_{h}$ is the horizontal gradient operator, $h_{s,t}$ is the vertical bottom velocity representing displacement of the impermeable substrate and $g$ is gravitational acceleration. Matching conditions at the water–sediment interface $z=-h+{\it\eta}_{2}$ consist of continuity of pressure and kinematic constraints for each layer. After linearizing with respect to the pressure perturbation and small interface displacement, the resulting conditions are

(3.4) $$\begin{eqnarray}\left\{\begin{array}{@{}l@{}}(R-1)g{\it\eta}_{2}={\it\Phi}_{t}-R\mathscr{Q}_{t}\quad \\ W_{w}=W_{s}=(-h+{\it\eta}_{2})_{t}\quad \end{array}\right.\text{at}\;z=-h\end{eqnarray}$$

where $R={\it\rho}_{s}/{\it\rho}$ . The normal velocities at the interface inside the water column, $W_{w}$ , and sediment, $W_{s}$ , are given by:

(3.5) $$\begin{eqnarray}\left\{\begin{array}{@{}l@{}}W_{w}={\it\Phi}_{z}+\boldsymbol{{\rm\nabla}}_{h}h\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{h}{\it\Phi}\quad \\ W_{s}=\mathscr{Q}_{z}+\boldsymbol{{\rm\nabla}}_{h}h\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{h}\mathscr{Q}\quad \end{array}\right.\text{at}\;z=-h.\end{eqnarray}$$

3.2. Derivation of the mild-slope equation

A non-dimensionalization of the sediment layer equation yields a parameter ${\it\epsilon}={\it\omega}{\it\nu}_{s}/c_{s}^{2}$ characterizing the size of the damping term relative to the undamped wave equation, where ${\it\omega}$ represents angular wave frequency. For the cases considered here, ${\it\epsilon}=O(10^{-1})$ , and we treat the damping effect as a perturbation to the leading-order inviscid problem. The mild-slope equation is developed using the eigenfunctions for the two-layer, inviscid, compressible fluid problem with constant layer thicknesses $h$ and $a=h_{s}-h$ and a free surface. Retaining damping in the lower layer eigenfunctions, as in Silva et al. (Reference Silva, Salles and Govaere2003), would eliminate the leading-order damping term in the model but would involve complex-valued separation constants and resulting complex-valued model coefficients.

The upper and lower layer potentials may be expanded according to

(3.6) $$\begin{eqnarray}{\it\Phi}\left(x,y,z,t\right)=\mathop{\sum }_{n=0}^{\infty }{\it\Phi}_{n}\left(x,y,z,t\right)=\mathop{\sum }_{n=0}^{\infty }{\it\psi}_{n}\left(x,y,t\right)M_{n}(z)\end{eqnarray}$$

for the water column and

(3.7) $$\begin{eqnarray}\mathscr{Q}\left(x,y,z,t\right)=\mathop{\sum }_{n=0}^{\infty }\mathscr{Q}_{n}\left(x,y,z,t\right)=\mathop{\sum }_{n=0}^{\infty }{\it\psi}_{n}\left(x,y,t\right)N_{n}(z)\end{eqnarray}$$

in the sediment layer. For simple harmonic motion with frequency ${\it\omega}$ and separation constant ${\it\beta}^{2}$ in the vertical, the eigenfunctions $M_{n}(z)$ and $N_{n}(z)$ for the upper and lower layers are given by

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle M_{n}=\frac{(1-{\it\lambda}_{n}T_{n})\cosh ({\it\beta}_{w,n}(h+z))+({\it\lambda}_{n}-T_{n})\sinh ({\it\beta}_{w,n}(h+z))}{(1-{\it\lambda}_{n}T_{n})\cosh ({\it\beta}_{w,n}h)+({\it\lambda}_{n}-T_{n})\sinh ({\it\beta}_{w,n}h)}, & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle N_{n}=\frac{({\it\lambda}_{n}-T_{n})\cosh {\it\beta}_{s,n}(h_{s}+z)}{{\it\alpha}_{n}\sinh ({\it\beta}_{s,n}a)[(1-{\it\lambda}_{n}T_{n})\cosh ({\it\beta}_{w,n}h)+({\it\lambda}_{n}-T_{n})\sinh ({\it\beta}_{w,n}h)]}, & \displaystyle\end{eqnarray}$$

where $T_{n}=\tanh ({\it\beta}_{w,n}h)$ , ${\it\lambda}_{n}={\it\omega}^{2}/g{\it\beta}_{w,n}$ and ${\it\alpha}_{n}={\it\beta}_{s,n}/{\it\beta}_{w,n}$ . The separation constants ${\it\beta}_{w,n}$ and ${\it\beta}_{s,n}$ for water and sediment layers respectively are given by

(3.10a,b ) $$\begin{eqnarray}{\it\beta}_{w,n}^{2}=k_{n}^{2}-\frac{{\it\omega}^{2}}{c^{2}};\quad {\it\beta}_{s,n}^{2}=k_{n}^{2}-\frac{{\it\omega}^{2}}{c_{s}^{2}}\end{eqnarray}$$

where $k_{n}$ is the wavenumber. The eigenfunctions $M_{n}$ and $N_{n}$ form a complete Sturm–Liouville basis subject to the orthogonality constraint $I_{mn}+RK_{mn}=0;m\neq n$ , where

(3.11a,b ) $$\begin{eqnarray}I_{mn}=\int _{-h}^{0}M_{m}(z)M_{n}(z)\text{d}z;\quad K_{mn}=\int _{-h_{s}}^{-h}N_{m}(z)N_{n}(z)\text{d}z.\end{eqnarray}$$

$M_{n}$ and $N_{n}$ take the following values at the vertical boundaries:

(3.12) $$\begin{eqnarray}\left\{\begin{array}{@{}c@{}}M_{n}=1,\quad M_{n}^{\prime }={\it\omega}^{2}/g\quad \text{at}~z=0\\ M_{n}^{\prime }=N_{n}^{\prime }\quad \text{at}~z=-h\\ N_{n}^{\prime }=0\quad \text{at}~z=-h_{s}.\end{array}\right.\end{eqnarray}$$

The dispersion relation governing ${\it\beta}_{w,n}$ and ${\it\beta}_{s,n}$ is given by

(3.13) $$\begin{eqnarray}{\it\lambda}_{n}^{2}(R+{\it\alpha}_{n}T_{n}\hat{T}_{n})-{\it\lambda}_{n}R(T_{n}+{\it\alpha}_{n}\hat{T}_{n})+(R-1){\it\alpha}_{n}T_{n}\hat{T}_{n}=0,\end{eqnarray}$$

where $\hat{T}_{n}=\tanh ({\it\beta}_{s,n}a)$ . Equation (3.13) is a quartic system in ${\it\omega}$ describing a doubly infinite set of surface waves (with horizontal displacements in phase at the layer interface) and internal waves (with horizontal displacements 180° out of phase). The real roots of the dispersion relation ( $n=0$ ) are responsible for the primary surface and internal gravity waves, while the imaginary separation variables for $n\geqslant 1$ describe both progressive and spatially decaying hydro-acoustic modes. Due to the presence of damping in the real problem, the internal or interfacial wave modes are expected to be rapidly damped. If $a=0$ , then there is no sediment layer and (3.13) reduces to the classical hydro-acoustic dispersion relation given by ${\it\lambda}_{n}=\tanh ({\it\beta}_{w,n}h)$ .

We multiply the governing equation for each layer by each member of the set of eigenfunctions and integrate over the layer depth, giving

(3.14a,b ) $$\begin{eqnarray}\text{I}^{m}=\int _{-h}^{0}M_{m}L({\it\Phi})\text{d}z=0;\quad \text{II}^{m}=\int _{-h_{s}}^{-h}N_{m}L_{s}(\mathscr{Q})\text{d}z=0.\end{eqnarray}$$

Each expression is manipulated by introducing the expansions (3.6) and (3.7) and using the Leibniz rule and appropriate boundary conditions. Subsequently, we neglect second-order terms in the interfacial and substrate slope, staying within the classic mild-slope framework. The expressions for the two layers become

(3.15) $$\begin{eqnarray}\displaystyle \text{I}^{m} & = & \displaystyle \mathop{\sum }_{n=0}^{\infty }\left\{\boldsymbol{{\rm\nabla}}_{h}\boldsymbol{\cdot }[I_{mn}\boldsymbol{{\rm\nabla}}_{h}{\it\psi}_{n}]-\left(\left[\frac{I_{mn}}{c^{2}}+\frac{1}{g}\right]{\it\psi}_{n,t}\right)_{,t}-J_{mn}{\it\psi}_{n}\right.\nonumber\\ \displaystyle & & \displaystyle -\,M_{m}(-h)[M_{n}^{\prime }(-h){\it\psi}_{n}+M_{n}(-h)\boldsymbol{{\rm\nabla}}_{h}h\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{h}{\it\psi}_{n}]\nonumber\\ \displaystyle & & \displaystyle \left.+\,\frac{1}{c^{2}}M_{m}(-h)M_{n}(-h)h_{t}{\it\psi}_{n,t}\right\}\end{eqnarray}$$

for the water layer and

(3.16) $$\begin{eqnarray}\displaystyle \text{II}^{m} & = & \displaystyle \mathop{\sum }_{n=0}^{\infty }\left\{\boldsymbol{{\rm\nabla}}_{h}\boldsymbol{\cdot }[K_{mn}\boldsymbol{{\rm\nabla}}_{h}{\it\psi}_{n}]-\left(\frac{K_{mn}}{c_{s}^{2}}{\it\psi}_{n,t}\right)_{,t}-L_{mn}{\it\psi}_{n}-2{\it\epsilon}\frac{{\it\omega}}{c_{s}^{2}}K_{mn}{\it\psi}_{n,t}\right.\nonumber\\ \displaystyle & & \displaystyle +\,N_{m}(-h)[N_{n}^{\prime }(-h){\it\psi}_{n}+N_{n}(-h)\boldsymbol{{\rm\nabla}}_{h}h\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{h}{\it\psi}_{n}]\nonumber\\ \displaystyle & & \displaystyle -\,N_{m}(-h_{s})[N_{n}^{\prime }(-h_{s}){\it\psi}_{n}+N_{n}(-h_{s})\boldsymbol{{\rm\nabla}}_{h}h\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{h}{\it\psi}_{n}]\nonumber\\ \displaystyle & & \displaystyle \left.-\,\frac{1}{c_{s}^{2}}\left[N_{m}(-h)N_{n}(-h)h_{t}-N_{m}(-h_{s})N_{n}(-h_{s})h_{s,t}\right]{\it\psi}_{n,t}\right\}\end{eqnarray}$$

for the sediment layer, where we have introduced the approximation ${\rm\nabla}^{2}\mathscr{Q}\approx (1/c_{s}^{2})\mathscr{Q}_{tt}\approx -({\it\omega}^{2}/c_{s}^{2})\mathscr{Q}$ in the damping term. Here, $J_{mn}$ and $L_{mn}$ in (3.15) and (3.16) are given by

(3.17) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle J_{mn}=\int _{-h}^{0}M_{m}^{\prime }(z)M_{n}^{\prime }(z)\text{d}z=\frac{{\it\omega}^{2}}{g}-{\it\beta}_{w,n}^{2}I_{mn}-M_{m}(-h)M_{n}^{\prime }(-h)\\ \displaystyle L_{mn}=\int _{-h_{s}}^{-h}N_{m}^{\prime }(z)N_{n}^{\prime }(z)\text{d}z=-{\it\beta}_{s,n}^{2}K_{mn}+N_{m}(-h)N_{n}^{\prime }(-h)-N_{m}(-h_{s})N_{n}^{\prime }(-h_{s}).\end{array}\right\}\end{eqnarray}$$

Combining (3.15) and (3.16) according to $\text{I}^{m}+R\text{II}^{m}=0$ , in order to take advantage of orthogonality within the spatial derivative terms, and making use the interfacial kinematic and dynamic boundary conditions, we obtain the desired mild-slope equation

(3.18) $$\begin{eqnarray}\displaystyle (I_{2}^{m}{\it\psi}_{m,t})_{t}-\boldsymbol{{\rm\nabla}}_{h}\boldsymbol{\cdot }[I_{1}^{m}\boldsymbol{{\rm\nabla}}_{h}{\it\psi}_{m}]+[{\it\omega}^{2}I_{2}^{m}-k_{m}^{2}I_{1}^{m}]{\it\psi}_{m}+2R{\it\epsilon}\frac{{\it\omega}}{c_{s}^{2}}K_{n}{\it\psi}_{m,t}=D_{1}^{m}h_{t}+D_{2}^{m}h_{s,t}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where we have further introduced the approximation ${\it\psi}_{n,tt}=-{\it\omega}^{2}{\it\psi}_{n};n\neq m$ in order to eliminate an apparent coupling of the individual model equations arising in the time derivatives. Model coefficients are given by

(3.19) $$\begin{eqnarray}\displaystyle & I_{1}^{m}=I_{mm}+RK_{mm}, & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle I_{2}^{m}=\frac{I_{mm}}{c^{2}}+R\frac{K_{mm}}{c_{s}^{2}}+\frac{1}{g}, & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & D_{1}^{m}=-[M_{m}(-h)-RN_{m}(-h)], & \displaystyle\end{eqnarray}$$

(3.22) $$\begin{eqnarray}\displaystyle & D_{2}^{m}=-RN_{m}(-h_{s}), & \displaystyle\end{eqnarray}$$

where $I_{mm}(x,y,t)$ and $K_{mm}(x,y,t)$ are given by

(3.23) $$\begin{eqnarray}\displaystyle & \displaystyle I_{mm}=\int _{-h}^{0}M_{m}^{2}\text{d}z=\frac{h}{2G_{1}^{m}}[(1-2{\it\lambda}_{m}T_{m}+{\it\lambda}_{m}^{2})+(1-{\it\lambda}_{m}^{2})G_{1}^{m}], & \displaystyle\end{eqnarray}$$

(3.24) $$\begin{eqnarray}\displaystyle & \displaystyle K_{mm}=\int _{-h_{s}}^{-h}N_{m}^{2}\text{d}z=\frac{h({\it\lambda}_{m}-T_{m})^{2}}{2{\it\alpha}_{m}^{3}T_{m}\hat{T}_{m}}\frac{1+G_{2}^{m}}{G_{1}^{m}}, & \displaystyle\end{eqnarray}$$

with $G_{1}^{m}=2{\it\beta}_{w,m}h/\!\sinh (2{\it\beta}_{w,m}h)$ and $G_{2}^{m}=2{\it\beta}_{s,m}a/\!\sinh (2{\it\beta}_{s,m}a)$ . Equation (3.18) is the final form of the hyperbolic mild-slope equation for weakly compressible fluid overlying a sediment viscous layer (MSEDWC). The time dependence of $I_{2}^{m}$ in (3.18) is formally required to obtain a correct energy balance in the time-dependent medium, but in practice is too weak to affect numerical results noticeably. The elliptic version for purely harmonic motion is obtained by taking the Fourier transform of (3.18) and is given by

(3.25) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}_{h}\boldsymbol{\cdot }[I_{1}^{m}\boldsymbol{{\rm\nabla}}_{h}{\it\Psi}_{m}]+\left[k_{m}^{2}I_{1}^{m}-2\text{i}R{\it\epsilon}\frac{{\it\omega}}{c_{s}^{2}}K_{m}\right]{\it\Psi}_{m}=-\text{i}{\it\omega}D_{1}^{m}H-\text{i}{\it\omega}D_{2}^{m}H_{s},\end{eqnarray}$$

where ${\it\psi}_{n}(x,y,t)={\it\Psi}_{n}(x,y)\text{e}^{\text{i}{\it\omega}t}$ , $h(x,y,t)=H(x,y)\text{e}^{\text{i}{\it\omega}t}$ and $h_{s}(x,y,t)=H_{s}(x,y)\text{e}^{\text{i}{\it\omega}t}$ . For specific conditions (3.18) and (3.25) give the same solution as that previously proposed by other authors; for example: if $R=1$ or $a=0$ , then there is no sediment layer, and (3.18) and (3.25) reduce to the equation given by Sammarco et al. (Reference Sammarco, Cecioni, Bellotti and Abdolali2013).

Figure 3. Results for the free-surface elevation time series at 96 km from the tsunamigenic source according to the 3D (light grey) and depth-integrated (black) models in a constant depth, $h_{p}=2200~\text{m}$ , $a=1000~\text{m}$ , $c=1500~\text{m}~\text{s}^{-1}$ , $c_{s}=2000~\text{m}~\text{s}^{-1}$ , ${\it\rho}=1000~\text{kg}~\text{m}^{-3}$ , ${\it\rho}_{s}=1850~\text{kg}~\text{m}^{-3}$ , for a unit source area with semi-length $b=112~\text{km}$ and rise time ${\it\tau}=1~\text{s}$ . (a,c) Time series and (b,d) corresponding spectra with ${\it\gamma}_{1}={\it\omega}_{1}/2{\rm\pi}=0.15~\text{Hz}$ . (a,b) ${\it\mu}_{s}=0$ ; (c,d) ${\it\mu}_{s}=2\times 10^{8}~\text{Pa}~\text{s}$ .

Figure 4. The case of varying sea bottom and sediment thickness with a tsunamigenic source at the shallower part. (a) The computational domain. (b,d) Time series of free-surface elevation and (c,e) corresponding spectra at 400 km from the tsunamigenic source from the 3D (light grey) and depth-integrated models (black), $c=1500~\text{m}~\text{s}^{-1}$ , $c_{s}=2000~\text{m}~\text{s}^{-1}$ , ${\it\rho}=1000~\text{kg}~\text{m}^{-3}$ , ${\it\rho}_{s}=1850~\text{kg}~\text{m}^{-3}$ and ${\it\mu}_{s}=0$ for a unit source area with semi-length $b=15~\text{km}$ and rise time ${\it\tau}=1~\text{s}$ . (b,c) Impermeable sea bottom; (d,e) coupled model.

4. Sample computations

Sample computations have been carried out to verify whether the model (3.18) can be safely applied in place of a more computationally expensive 3D treatment. We present results for three different domains consisting of vertical sections in $x,z$ through laterally uniform domains with no $y$ -dependence, the first with a constant water depth and sediment thickness and the other two with varying water depth and sediment thickness. Frequency bands of width ${\rm\Delta}f=0.02~\text{Hz}$ have been selected to discretize the forcing spectrum. For the first case, the numerical solvers are applied on a computational domain 200 km long; given the symmetry of the problem about the middle of the earthquake ( $x=0$ ), computations are undertaken only for half of the physical domain. The Sommerfeld radiation condition is applied at the open end of the domain, so that the waves leave the domain freely. At $x=0$ , a fully reflective boundary condition is used in order to preserve symmetry. To correctly reproduce the wavefield, the maximum mesh size is 200 m, for a total of 1000 elements in the case of the depth-integrated model (3.18) and more than 30 000 elements for the 3D model (3.1) and (3.2). The time step is $t=0.1~\text{s}$ and the computational time to reproduce 500 s of real-time simulation was approximately 10 min for (3.18) and approximately 3 h for (3.1) and (3.2); a computer equipped with an i7 3.2 GHz CPU and 64 GB RAM has been used. The simplified earthquake effect is modelled as a displacement in the vertical direction of the bottom with duration ${\it\tau}=1~\text{s}$ and bottom velocity defined by (2.3); the results are presented in figure 3 in terms of the free-surface elevation ${\it\eta}$ and the corresponding spectrum $\tilde{{\it\eta}}$ . Results are shown for a virtual surface gauge at $x=96~\text{km}$ over a 112 km semi-fault where water depth and sediment thickness are $h_{p}=2200$ and $a=1000~\text{m}$ respectively. In 3(a,b) the results relate to the case of inviscid sediment. The 3D model (light grey) and depth-integrated model (black) are in optimal agreement. The peak frequency is ${\it\gamma}=0.15~\text{Hz}$ , corresponding to the first cutoff frequency for the coupled system defined by (2.2). The hydro-acoustic waves remain at the same order as the generated waves until the end of computations. In 3(c,d), we add the dissipation term into the equation ( ${\it\mu}_{s}=2\times 10^{8}~\text{Pa}~\text{s}$ ). The generated hydro-acoustic waves are attenuated gradually during 500 s.

Figure 5. The case of varying sea bottom and sediment thickness with a tsunamigenic source at the deeper part. (a) The computational domain. (b,d) Time series of free-surface elevation and (c,e) corresponding spectra at point $A$ , 100 km from the tsunamigenic source at 3.5 km water depth (black), and point $B$ , 400 km from the tsunamigenic source at 2 km water depth (light grey), obtained from depth-integrated models. The water and sediment characteristics are the same as in figure 4. (b,c) Impermeable sea bottom; (d,e) coupled model.

In the second case, a varying sea water depth overlying a varying sedimentary layer is considered; the domain’s geometry, depicted in the upper plot of figure 4, has a region of length 220 km with a depth of 2 km in each layer, a region of length 50 km with a sloping bottom, and another region of length 230 km with a constant water depth of 3.5 km and a sediment thickness of 0.5 km. Over the entire 500 km domain length, the total depth $h_{s}=4~\text{km}$ . Ground motion occurs in a patch of length 15 km at the left edge of the shallow area, as shown in figure 4(a); it moves vertically with a bottom velocity and total displacement of 1 m given by (2.3). The maximum mesh size is again 200 m, for a total of 2500 elements in the case of the depth-integrated model (3.18) and 130 000 triangular elements for the 3D model ((3.1) and (3.2)). The time step and discretization of the spectra are the same as in the constant-depth case. The computational time to reproduce $1000~\text{s}$ of real-time simulation was approximately 20 min for (3.18) and approximately 5 h for (3.1) and (3.2), using the same computer as in the previous simulation. The results are presented in figure 4 in terms of time series of the free-surface elevation ${\it\eta}$ and the corresponding spectra $\tilde{{\it\eta}}$ , at a distance $x=400~\text{km}$ from the moving sea bed area. The two time series are in good agreement, both in terms of amplitude and modulation of the signal. The comparison results show that the peak frequency shifts from $f_{1}=0.1875~\text{Hz}$ for an impermeable bottom to ${\it\gamma}_{1}=0.138~\text{Hz}$ for the coupled water and sediment model.

In the third case, model performance has been investigated for the case of an earthquake on the deeper part of continental shelf in order to reveal the transmission and reflection properties of the hydro-acoustic wavefield in deep and shallow waters. The domain is depicted in figure 5(a), where the earthquake occurs in the deeper area (3.5 km water depth over a 0.5 km sediment) with the same physical characteristics for the water, sediment and tsunamigenic source as in the previous cases. The results are presented in figure 5 in terms of time series of the free-surface elevation ${\it\eta}$ and the corresponding spectra $\tilde{{\it\eta}}$ at point $A$ , 100 km from the epicentre at 3.5 km water depth (black), and point $B$ , 400 km from the tsunamigenic source at 2 km water depth (light grey), obtained from depth-integrated models for an impermeable (b,c) and a permeable (d,e) sea bottom. It can be seen from the model results that the hydro-acoustic waves cannot propagate upslope. Hydro-acoustic wave frequencies lower than the corresponding cutoff frequency of the observatory depth have been filtered. As a result, the wave amplitudes decreased from point $A$ to $B$ . The waves reflected from the slope are superimposed on the arriving wave trains and change the wave packet modulation. This result indicates the importance of deep sea observatories for hydro-acoustic wave detection (Abdolali et al. Reference Abdolali, Cecioni, Bellotti and Sammarco2014).