2.1. Floating plate in shallow water

We first consider an infinite elastic flexible plate floating on surface of the fluid (figure 1). Assuming that no separation occurs, surface elevation is the same as displacement of the plate. The kinematic and dynamic surface boundary conditions resulting from the Kirchhoff–Love plate theory are given as

(2.11)$$\begin{gather} \frac{{\partial} \eta^*}{{\partial} t^*} + (\boldsymbol{u}^* \boldsymbol{\cdot} \boldsymbol{\nabla}^*)\eta^* = w^*,\quad {\rm at}\ z^*=\eta^*(x^*, y^*, t^*), \end{gather}$$

(2.12)$$\begin{gather}\alpha^*\,\frac{{\partial}^2 \eta^*}{{\partial} {t^*}^2} + {\nabla^*}^2(\beta^*\,{\nabla^*}^2\eta^*) + P^*_a - P^* = 0,\quad {\rm at}\ z^*=\eta^*(x^*, y^*, t^*), \end{gather}$$

in which ${\nabla ^*}^2={ {\partial }^2}/{ {\partial } {x^*}^2}+{ {\partial }^2}/{ {\partial } {y^*}^2}$ is the Laplacian operator, $P^*_a$ is the ambient air pressure on top of the plate, assumed to be zero, $\alpha ^*(x^*,y^*)=\rho _{plt}\, d^*(x^*,y^*)$ is the mass of the plate per unit area, $\beta ^*(x^*,y^*)={E{d^*}^3}/{12(1-\nu ^2)}$ is the flexural stiffness of the plate, $d^*$ is the thickness of the plate, assumed to be negligible, and $\rho _{plt}$, $E$ and $\nu$ are plate density, Young's modulus and Poisson's ratio, respectively. Introducing the additional non-dimensional variables

(2.13a,b)\begin{equation} \alpha=\alpha^*/(\rho H_0),\quad \beta=\beta^*/(\rho g{L_0}^4), \end{equation}

(2.11) and (2.12) can be non-dimensionalized as

(2.14)$$\begin{gather} \frac{{\partial} \eta}{{\partial} t} + \epsilon(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})\eta -w = 0,\quad {\rm at}\ z=\epsilon\,\eta(x, y, t), \end{gather}$$

(2.15)$$\begin{gather}\mu^2\alpha\,\frac{{\partial}^2 \eta}{{\partial} {t}^2} + {\nabla}^2(\beta\,{\nabla}^2\eta) -P = 0,\quad {\rm at}\ z=\epsilon\,\eta(x, y, t). \end{gather}$$

An expansion method is employed to expand non-dimensional variables with respect to $\epsilon$ for the derivation of Boussinesq-type equations, which is similar to Peregrine's work (Peregrine Reference Peregrine1967) but with major modifications to include the elastic and inertial effects of the plate. Details of the mathematical derivation are given in Appendix A. The resulting governing equations are presented in the following subsections.

2.1.1. The zeroth-order solutions and first-order equations

The zeroth-order solutions result from the still-water condition with no displacement:

(2.16a–c)\begin{equation} \boldsymbol{u}_0=\boldsymbol{Q}_0=\boldsymbol{0},\quad w_0=0\quad {\rm and}\quad P_0=-z/\epsilon, \end{equation}

in which $\boldsymbol {Q}$ is the horizontal flux of the flow, and the subscript indicates the associated order with respect to $\epsilon$ in the expansion series.

The first-order linear governing equation that describes FG waves propagating in a floating plate with variable stiffness over a varying bathymetry are derived as

(2.17)\begin{equation} \frac{{\partial}^2 \eta_1}{{\partial} t^2} - \boldsymbol{\nabla} \boldsymbol{\cdot}(H\,\boldsymbol{\nabla}\eta_1) - \boldsymbol{\nabla} \boldsymbol{\cdot}\{H\,\boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\eta_1)]\} = 0. \end{equation}

Equation (2.17) is applicable to mild variations of water depth and plate stiffness, and can be regarded as the ‘mild-variation’ equation of FG waves, which is similar to the time-dependent mild-slope equation of shallow-water waves in Bin (Reference Bin1994). Equation (2.17) reduces to the standard mild-slope equation when $\beta =0$.

Once (2.17) is solved for $\eta _1$, other first-order variables can be obtained: $P_1$ can be solved from (A21); $\boldsymbol {u}_1$ and $\boldsymbol {Q}_1$ can be computed using (A22) and (A24); and $w_1$ can be calculated using the following equation that results from the integration of (A15),

(2.18)\begin{equation} w_1 =-z\,\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_1 - \boldsymbol{\nabla} \boldsymbol{\cdot}(H\boldsymbol{u}_1), \end{equation}

which satisfies the boundary conditions of (A18) and (A20).

2.1.2. Boussinesq-type equations

The higher-order nonlinear governing equations are derived in Appendix A, and the resulting governing equations are presented as

(2.19)\begin{gather} \begin{aligned} & \frac{{\partial} \bar{\boldsymbol{u}}}{{\partial} t} + \epsilon(\bar{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}) \bar{\boldsymbol{u}} + \boldsymbol{\nabla}\eta + \boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\eta)] + \mu^2\,\boldsymbol{\nabla}\left(\alpha\,\frac{{\partial}^2\eta}{{\partial} t^2}\right) \nonumber\\ &\quad = \mu^2\,\frac{H}{2}\,\frac{{\partial}}{{\partial} t} \boldsymbol{\nabla}[\boldsymbol{\nabla}\boldsymbol{\cdot}(H\bar{\boldsymbol{u}})] - \mu^2\,\frac{H^2}{6}\,\frac{{\partial}}{{\partial} t}\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{u}}), \end{aligned} \end{gather}

(2.20)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}[(H + \epsilon\eta)\bar{\boldsymbol{u}}] + \frac{{\partial} \eta}{{\partial} t}=0, \end{gather}

in which $\bar {\boldsymbol {u}}$ is the mean horizontal velocity defined by (A39).

Equations (2.19) and (2.20) are the Boussinesq-type equations for FG waves of a floating plate in shallow water, in which the elastic and inertial effects of the plate are represented respectively by the terms $\boldsymbol {\nabla }[\nabla ^2(\beta \,\nabla ^2\eta )]$ and $\mu ^2\,\boldsymbol {\nabla }(\alpha ({ {\partial }^2\eta }/{ {\partial } t^2}))$. These two coupled equations are used to solve for $\bar {\boldsymbol {u}}(x,y,t)$ and $\eta (x,y,t)$. They are reduced to the classic Boussinesq equations for shallow-water waves in Peregrine (Reference Peregrine1967) when $\alpha =\beta =0$. The Boussinesq-type equations can be represented by (A43) and (A44) in a dimensional form as well.

2.2. Submerged plate in shallow water

Consider an infinite elastic flexible plate submerged in the fluid (figure 2). The free-surface boundary conditions are applicable, whose non-dimensional forms are presented as

(2.21)$$\begin{gather} \frac{{\partial} \hat{\eta}}{{\partial} t} + \epsilon(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \hat{\eta} - w = 0,\quad {\rm at}\ z=\epsilon\,\hat{\eta}(x, y, t), \end{gather}$$

(2.22)$$\begin{gather}P = P_a = 0,\quad {\rm at}\ z=\epsilon\,\hat{\eta}(x, y, t), \end{gather}$$

where $\hat {\eta }=\hat {\eta }^*/A_0$ denotes non-dimensional free-surface elevation, which is different from the floating plate displacement $\eta$ as they satisfy different dynamic boundary conditions.

Figure 2. Submerged elastic plate over varying bathymetry.

Due to the presence of the submerged plate, an additional set of boundary conditions is introduced to represent the kinematics and dynamics of the plate:

(2.23)$$\begin{gather} w^+- (\boldsymbol{u}^+ \boldsymbol{\cdot} \boldsymbol{\nabla})(\epsilon\zeta-H_I) = w^\unicode{x2013}- (\boldsymbol{u}^- \boldsymbol{\cdot}\boldsymbol{\nabla})(\epsilon\zeta-H_I) = \frac{{\partial} \zeta}{{\partial} t},\quad {\rm at}\ z=\epsilon\zeta-H_I, \end{gather}$$

(2.24)$$\begin{gather}\mu^2\alpha\,\frac{{\partial}^2 \zeta}{{\partial} {t}^2} + {\nabla}^2(\beta\,{\nabla}^2\zeta) + P^+- P^-= 0,\quad {\rm at}\ z=\epsilon\zeta-H_I, \end{gather}$$

in which $H_I=H_I^*/H_0$ is non-dimensional submerged depth, $w^+ = w|_{z=-H_I+\epsilon \zeta ^+}$, $w^- = w|_{z=-H_I+\epsilon \zeta ^-}$, $\boldsymbol {u}^+ = \boldsymbol {u}|_{z=-H_I+\epsilon \zeta ^+}$ and $\boldsymbol {u}^- = \boldsymbol {u}|_{z=-H_I+\epsilon \zeta ^-}$ are fluid velocities at the upper and lower surfaces of the plate, and $P^+ = P|_{z=-H_I+\epsilon \zeta ^+}$ and $P^- = P|_{z=-H_I+\epsilon \zeta ^-}$ are fluid pressures on the upper and lower plate surfaces.

The flow field is separated by the submerged plate into two layers, each of which has its own boundary conditions but conforms to the same equations of motion as (2.6)–(2.9). We define the horizontal flux of the upper flow as $\boldsymbol {Q}_I=\int _{-H_I+\epsilon \zeta ^+}^{\epsilon \hat {\eta }} \boldsymbol {u}\,{\rm d} z$, and the lower horizontal flux as $\boldsymbol {Q}_{II}=\int _{-H}^{-H_I+\epsilon \zeta ^-} \boldsymbol {u}\,{\rm d} z$. Therefore, integrating (2.8) from $z=-H_I+\epsilon \zeta ^+$ to $z=\epsilon \hat {\eta }$, and from $z=-H$ to $z=-H_I+\epsilon \zeta ^-$, respectively, we have

(2.25)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{Q}_I + \frac{{\partial} \hat{\eta}}{{\partial} t} - \frac{{\partial} \zeta}{{\partial} t} = 0, \end{gather}$$

(2.26)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{Q}_{II} + \frac{{\partial} \zeta}{{\partial} t} = 0, \end{gather}$$

in which the boundary conditions of $w$ from (2.10), (2.21) and (2.23) have been applied.

2.2.1. The zeroth-order solutions and first-order equations

The free surface elevation $\hat {\eta }$, the submerged plate displacement $\zeta$, and the horizontal fluxes $\boldsymbol {Q}_I$ and $\boldsymbol {Q}_{II}$ are similarly expanded with respect to $\epsilon$ as in Appendix A. The system of governing equations is then obtained by substituting the expanded variables into the equations of motion and boundary conditions, including (2.6), (2.7), (2.9), (2.10) and (2.21)–(2.26).

The zeroth-order solution of these governing equations is the still-water condition:

(2.27a–d)\begin{align} \boldsymbol{u}_0=\boldsymbol{Q}_{I0}=\boldsymbol{Q}_{II0}=\boldsymbol{0},\quad w_0=0,\quad P_0=-z/\epsilon\quad {\rm and}\quad P_0^+=P_0^-=-H_I/\epsilon. \end{align}

Following steps similar to those in Appendix A, we can derive the first-order linear equations describing propagation of the coupled surface gravity wave and SH wave in shallow water. The resulting equations are

(2.28)$$\begin{gather} \frac{{\partial}^2 \hat{\eta}_1}{{\partial} t^2} - \frac{{\partial}^2 \zeta_1}{{\partial} t^2} - \boldsymbol{\nabla} \boldsymbol{\cdot}(H_I\,\boldsymbol{\nabla}\hat{\eta}_1) = 0, \end{gather}$$

(2.29)$$\begin{gather}\frac{{\partial}^2 \zeta_1}{{\partial} t^2} - \boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\,\boldsymbol{\nabla}\hat{\eta}_1) - \boldsymbol{\nabla}\boldsymbol{\cdot}\{H_{II}\,\boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\zeta_1)]\}= 0, \end{gather}$$

in which $H_{II}=H-H_I$ is the non-dimensional distance from the submerged plate to the bottom. Equation (2.29) reduces to (2.17) when $H_I$ goes to zero and $\zeta _1$ becomes $\eta _1$. Similar to (2.17), then (2.28) and (2.29) are the ‘mild-variation’ equations of SH waves applicable to mild variations of $H_I$, $H_{II}$ and $\beta$. The first-order surface elevation and submerged plate displacement are solved using (2.28) and (2.29), and other first-order variables can be computed subsequently.

2.2.2. Boussinesq-type equations

The higher-order equations can be derived similarly to the derivation in Appendix A. Because the submerged plate separates the flow into two layers, the number of unknown variables is doubled. We define the mean horizontal velocities of the upper and lower flows as

(2.30a,b)\begin{equation} \bar{\boldsymbol{u}}_I = \frac{\boldsymbol{Q}_I}{H_I+\epsilon\hat{\eta}-\epsilon\zeta}\quad {\rm and}\quad \bar{\boldsymbol{u}}_{II} = \frac{\boldsymbol{Q}_{II}}{H_{II} + \epsilon\zeta}. \end{equation}

The Boussinesq-type equations for the gravity–hydroelastic coupled wave system that consists of surface gravity waves and SH waves are derived in terms of mean horizontal velocities and surface displacements:

(2.31)\begin{gather} \frac{{\partial} \bar{\boldsymbol{u}}_I}{{\partial} t} + \epsilon(\bar{\boldsymbol{u}}_I\boldsymbol{\cdot}\boldsymbol{\nabla}) \bar{\boldsymbol{u}}_I + \boldsymbol{\nabla}\hat{\eta} = \mu^2\,\frac{H_I}{2}\,\frac{{\partial}}{{\partial} t} \boldsymbol{\nabla}[\boldsymbol{\nabla}\boldsymbol{\cdot}(H\bar{\boldsymbol{u}}_I)] - \mu^2\,\frac{H_I^2}{6}\,\frac{{\partial}}{{\partial} t} \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{u}}_I), \end{gather}

(2.32)\begin{gather} \begin{aligned} & \frac{{\partial} \bar{\boldsymbol{u}}_{II}}{{\partial} t} + \epsilon(\bar{\boldsymbol{u}}_{II}\boldsymbol{\cdot}\boldsymbol{\nabla})\bar{\boldsymbol{u}}_{II} + \boldsymbol{\nabla}\hat{\eta} + \boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\zeta)] + \mu^2\,\boldsymbol{\nabla}\left(\alpha\,\frac{{\partial}^2\zeta}{{\partial} t^2}\right) \nonumber\\ &\quad =\mu^2\,\frac{(H^2-H_I^2)}{2H_{II}}\,\frac{{\partial}}{{\partial} t} \boldsymbol{\nabla}[\boldsymbol{\nabla}\boldsymbol{\cdot}(H\bar{\boldsymbol{u}}_{II})] - \mu^2\, \frac{(H^3-H_I^3)}{6H_{II}}\,\frac{{\partial}}{{\partial} t}\boldsymbol{\nabla} (\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{u}}_{II}), \end{aligned}\end{gather}

(2.33)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}[(H_I+\epsilon\hat{\eta}-\epsilon\zeta)\bar{\boldsymbol{u}}_I] + \frac{{\partial} \hat{\eta}}{{\partial} t} - \frac{{\partial} \zeta}{{\partial} t} = 0, \end{gather}$$

(2.34)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}[(H_{II}+\epsilon\zeta)\bar{\boldsymbol{u}}_{II}] + \frac{{\partial} \zeta}{{\partial} t} = 0. \end{gather}$$

3.1. First-order analytical solution

Analytical solutions are highly desired because they enable straightforward computation and provide meaningful insights into the fluid–structure coupled dynamics. However, it is difficult to analytically solve the nonlinear Boussinesq-type equations. Instead, we derive the analytical solution to the first-order, linear governing equation, while numerically solving the nonlinear Boussinesq-type equations of FG waves.

The first-order, linear, non-dimensional governing equation of shallow-water FG waves is given by (2.17). Assuming that water depth $H$ and plate stiffness $\beta$ are constants that do not vary in space, (2.17) is reduced to

(3.1)\begin{equation} \frac{{\partial}^2 \eta_1}{{\partial} t^2} - H\,\frac{{\partial}^2 \eta_1}{{\partial} x^2} - H\beta\,\frac{{\partial}^6 \eta_1}{{\partial} x^6} = 0, \end{equation}

in which $x$ is the direction of wave propagation. Assume that the solution to (3.1) is periodic: $\eta _1=\eta _a \exp ({{\rm i}(kx\pm \sigma t)})$, in which $k$ and $\sigma$ are non-dimensional wavenumber and angular frequency, respectively, and $\eta _a$ is the amplitude determined by boundary and initial conditions. Substituting the periodic solution into (3.1), we can derive the first-order dispersion relation for FG waves in shallow water:

(3.2)\begin{equation} \sigma^2=Hk^2+H\beta k^6. \end{equation}

Equation (3.2) reduces to the dispersion relation of linear shallow-water waves when flexural stiffness $\beta$ goes to zero. Substituting $\gamma =k^2$ into (3.2), and solving the resulting cubic equation for the real root of $\gamma$ using Cardano's formula, we can compute the wavenumber for the given wave frequency, water depth and plate stiffness:

(3.3)\begin{equation} {k}^2=\sqrt[3]{\frac{{\sigma}^2}{2H\beta}-\sqrt{\left(\frac{{\sigma}^2}{2H\beta}\right)^2+ \left(\frac{1}{3\beta}\right)^3}}+\sqrt[3]{\frac{{\sigma}^2}{2H\beta}+ \sqrt{{\left(\frac{{\sigma}^2}{2H\beta}\right)}^2+{\left(\frac{1}{3\beta}\right)}^3}}. \end{equation}

Additionally, assume that the periodic horizontal velocity is $u_1=u_a \exp ({{\rm i}(kx\pm \sigma t)})$, where $u_a$ is the amplitude of the velocity. Substituting the periodic solutions of $\eta _1$ and $u_1$ into (A16), (A22) and (A24), and combining the resulting equations, we have

(3.4)\begin{equation} u_a=\frac{k+\beta k^5}{\sigma}\,\eta_a=\frac{\sigma}{Hk}\,\eta_a. \end{equation}

Equations (3.3) and (3.4) offer the complete periodic solutions to the first-order linear system of FG waves in shallow water.

3.2. Linear dispersion relations

Equation (3.2) presents the dispersion relation derived from the first-order equation, in which all the terms associated with $\epsilon$ and $\mu ^2$ in the Boussinesq-type equations are neglected. Now we derive the linear dispersion relation from the reduced equations of (2.19) and (2.20), in which only the nonlinear terms associated with $\epsilon$ are neglected. Similarly assuming constant depth and uniform plate, the reduced equations are

(3.5)$$\begin{gather} \frac{{\partial} \bar{u}_x}{{\partial} t} + \frac{{\partial} \eta}{{\partial} x} + \beta\,\frac{{\partial}^5 \eta}{{\partial} x^5} + \mu^2\alpha\,\frac{{\partial}^3\eta}{{\partial} x\,{\partial} t^2} - \mu^2\,\frac{H^2}{3}\,\frac{{\partial}^3 \bar{u}_x}{{\partial} x^2\,{\partial} t} = 0, \end{gather}$$

(3.6)$$\begin{gather}H\,\frac{{\partial} \bar{u}_x}{{\partial} x} + \frac{{\partial} \eta}{{\partial} t} = 0, \end{gather}$$

in which $\bar {u}_x$ is the $x$-component of $\bar {\boldsymbol {u}}$ in the propagating direction.

The linear dispersion relation is obtained by substituting harmonically periodic $\bar {u}_x$ and $\eta$ into (3.5) and (3.6):

(3.7)\begin{equation} \sigma^2=\frac{Hk^2+H\beta k^6}{1+\mu^2H\alpha k^2+\tfrac{1}{3}\mu^2k^2H^2}. \end{equation}

Similar to § 3.1, the linear analytic solutions of $\bar {u}_x$ and $\eta$ can be obtained by solving the linear dispersion relation (3.7).

The linear dispersion relation of FG waves for finite water depth was derived by Fox & Squire (Reference Fox and Squire1994) and Sahoo, Yip & Chwang (Reference Sahoo, Yip and Chwang1994), and it can be non-dimensionalized as

(3.8)\begin{equation} \sigma^2=(Hk^2-\mu^2H\alpha k^2\sigma^2+H\beta k^6)\,\frac{\tanh(\mu kH)}{\mu kH}. \end{equation}

When $\mu$ is small such that the [$1/2$]-order Padé approximation of $\tanh (\mu kH)$ gives ${\mu k H}/{(1+(\mu kH)^2/3)}$, (3.8) is reduced to

(3.9)\begin{equation} \sigma^2=(Hk^2-\mu^2H\alpha k^2\sigma^2+H\beta k^6)\,\frac{1}{1+\tfrac{1}{3}\mu^2k^2H^2}, \end{equation}

which is equivalent to (3.7) if $\sigma ^2$ is moved to the left-hand side of the equation. Furthermore, (3.7) and(3.8) both converge to (3.2) if the higher-order terms associated with $\mu ^2$ are neglected.

The phase and group velocities of FG waves can be derived using different linear dispersion relations for shallow water and finite depth, respectively. For the first-order shallow-water model, its phase and group velocities are derived from (3.2) as

(3.10)$$\begin{gather} C_{p\_first}^2 = \frac{\sigma^2}{k^2}=H+H\beta k^4, \end{gather}$$

(3.11)$$\begin{gather}C_{g\_first} = \frac{{\partial} \sigma}{{\partial} k}=C_{p\_first}\,\frac{1+3\beta k^4}{1+\beta k^4}. \end{gather}$$

For the linear shallow-water model, its phase and group velocities are derived from (3.7):

(3.12)$$\begin{gather} C_{p\_linear}^2 = \frac{H+H\beta k^4}{1+\mu^2H\alpha k^2+\tfrac{1}{3}\mu^2k^2H^2}, \end{gather}$$

(3.13)$$\begin{gather}C_{g\_linear} = C_{p\_linear}\Big(\frac{2\beta k^4}{1+\beta k^4}+\frac{1}{1+\mu^2 H\alpha k^2+ \tfrac{1}{3}\mu^2 k^2H^2}\Big)\!. \end{gather}$$

For the linear finite depth model, its phase and group velocities are derived from (3.8):

(3.14)$$\begin{gather} C_{p\_finite}^2 = \frac{(H+H\beta k^4)\,\dfrac{\tanh(\mu kH)}{\mu kH}}{1+\mu^2H\alpha k^2\, \dfrac{\tanh(\mu kH)}{\mu kH}}, \end{gather}$$

(3.15)$$\begin{gather}C_{g\_finite} = C_{p\_finite}\Big[\frac{2\beta k^4}{1+\beta k^4}+\frac{\mu k H +\tanh(\mu k H)- \mu k H\tanh^2(\mu k H)}{2\tanh(\mu k H)+2\mu\alpha k\tanh^2(\mu k H)}\Big]\!. \end{gather}$$

To evaluate the effectiveness of the proposed shallow-water models for various relative depth $\mu$, we compute the phase and group velocities of FG waves propagating in sea-ice (table 1) using the first-order, linear shallow-water and linear finite-depth dispersion relations. The phase and group velocities of linear water waves (Airy waves) are included for comparison as well. Figure 3 shows the changes of phase and group velocities against $\mu$ with different ice thicknesses when a typical wave period $T^*=10$ s is considered. The phase and group velocities predicted by the first-order dispersion relation generally agree well with those of the finite-depth model up to $\mu =0.05$ (1.7 % error in phase velocity for $d^*=1.0$ m), while the linear shallow-water model can be applied up to $\mu =0.2$ (2.0 % error in phase velocity for $d^*=1.0$ m). It is noticed that the elastic effects are prominent in shallow water when we compare the water waves with FG waves: the predictions by Airy wave theory are closer to FG waves in finite depth for the larger $\mu$.

Table 1. Physical parameters of floating sea-ice.

Figure 3. Comparison of phase and group velocities computed using different wave models with various $\mu$ ($T^*=10$ s).

The effects of wave frequency are also investigated. Figure 4 shows the phase and group velocities of water waves on surface and FG waves in sea-ice for different wave periods. The wavelength of FG waves increases with wave period but its velocities do not, and the group velocity of FG waves is greater than the phase velocity for higher frequencies. These characteristics are completely different from water waves. For pure bending waves, the phase velocity is $\sqrt [4]{\sigma ^2\beta /\alpha }$, which decreases with wave period, and the group velocity is twice as large as the phase velocity. The FG waves result from the combining effects of gravitation and elastic flexion, and it is noticed that the elastic effects are prominent for higher frequencies.

Figure 4. Comparison of phase and group velocities computed for different wave periods: (a,b) $d^*=0$ (water wave); (c,d) $d^*=0.3$ m (FG wave); (e,f) $d^*=1.0$ m (FG wave).

3.3. Nonlinear numerical solutions

The nonlinear Boussinesq-type equations of FG waves in shallow water are solved numerically. Equations (2.19) and (2.20) are computed using an implicit, time-marching scheme:

(3.16)\begin{gather} \eta_{s+1^*} = \eta_{s} - \Delta t\,\boldsymbol{\nabla} \boldsymbol{\cdot} [(\epsilon\eta_{s}+H)\bar{\boldsymbol{u}}_{s}], \end{gather}

(3.17)\begin{gather} \begin{aligned}& \frac{\bar{\boldsymbol{u}}_{s+1}-\bar{\boldsymbol{u}}_{s}}{\Delta t} +\epsilon(\bar{\boldsymbol{u}}_s\boldsymbol{\cdot}\boldsymbol{\nabla})\bar{\boldsymbol{u}}_s -\frac{\mu^2}{2}\,H\,\boldsymbol{\nabla}\left[\boldsymbol{\nabla}\boldsymbol{\cdot}\left(H\,\frac{\bar{\boldsymbol{u}}_{s+1}- \bar{\boldsymbol{u}}_s}{\Delta t}\right)\right] \nonumber\\ &\qquad +\frac{\mu^2}{6}\,H^2\,\boldsymbol{\nabla} \left[\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\bar{\boldsymbol{u}}_{s+1}-\bar{\boldsymbol{u}}_s}{\Delta t}\right)\right] \nonumber\\ &\quad =-\boldsymbol{\nabla}\left(\frac{\eta_{s}+\eta_{s+1^*}}{2}\right) - \boldsymbol{\nabla}\left\{\nabla^2\left[\beta\,\nabla^2\left(\frac{\eta_{s}+\eta_{s+1^*}}{2}\right)\right]\right\} \nonumber\\ &\qquad -\mu^2\,\boldsymbol{\nabla}\left(\alpha\,\frac{\eta_{s+1^*}-2\eta_{s}+\eta_{s-1}}{\Delta t^2}\right)\!, \end{aligned}\end{gather}

(3.18)\begin{gather} \eta_{s+1} = \eta_{s} - \Delta t\,\boldsymbol{\nabla} \boldsymbol{\cdot} [(\epsilon\eta_{s+1^*}+H)\bar{\boldsymbol{u}}_{s+1}], \end{gather}

in which $\Delta t$ is the time step size, subscript $s$ indicates the $s$th time step, and $\eta _{s+1^*}$ is the provisional value of $\eta _{s+1}$. Equation (3.17) is solved implicitly to update $\bar {\boldsymbol {u}}_{s+1}$. Application of the provisional step allows conditionally stable computation of (2.19) and (2.20) using the central difference schemes for spatial derivatives.

To demonstrate the nonlinear effects, FG waves propagating in sea-ice in shallow water are computed using the first-order linear analytical solution and the proposed nonlinear numerical scheme, respectively. The parameters used in these computations are presented in tables 1 and 2, and thickness 0.1 m is employed. By changing water depth and the amplitude of boundary excitation, we compare the linear analytical solutions and the nonlinear Boussinesq solutions based on different values of $\epsilon$ and $\mu$ that generally conform to the Boussinesq approximations: $O(\epsilon ) \approx O(\mu ^2)\ll O(1)$. Figure 5 shows that the nonlinear FG waves have smaller wavelength than their first-order linear counterparts.

Table 2. Simulation cases.

Figure 5. Comparison between the linear and nonlinear solutions of FG waves with various $\epsilon$ and $\mu$. (a) $\epsilon = 0.0010,\ \mu = 0.0309$; (b) $\epsilon = 0.0051,\ \mu = 0.0713$; (c) $\epsilon = 0.0100,\ \mu = 0.1000$; (d) $\epsilon = 0.0204, \mu = 0.1428$.

The shoaling effect of FG waves is illustrated by the nonlinear numerical solutions. Figure 6 presents the snapshots of nonlinear FG waves in sea-ice with different thicknesses propagating over decreasing water depth. The wavelength reduces in all three cases during shoaling, regardless of the different ice thickness, as predicted by the linear dispersion relation. However, the amplitude of FG waves does not necessarily increase during shoaling. For flexural deflection, its elastic strain energy is associated with the curvature of deflection, which increases with decreasing wavelength if wave amplitude remains. When the elastic effect is dominant in an FG wave system (with relatively large plate stiffness), its amplitude must decrease with decreasing wavelength during shoaling for the conversation of energy. Recalling the relation between wave velocities and wave frequency shown in figure 4, this is another unique behaviour of shallow FG waves that is different from shallow-water waves, demonstrating the combining effects of gravitation and elastic flexion.

Figure 6. Shoaling effect of nonlinear FG waves. (a) $d^{\ast} = 0.1\ {\rm m}$; (b) $d^{\ast} = 0.3\ {\rm m}$; (c) $d^{\ast} = 0.5\ {\rm m}$; (d) Varying depth.

3.4. Solitary wave solution

A solitary wave solution is derived for nonlinear FG waves in shallow water. Similarly to the derivation of linear analytical solutions, we assume that water depth $H$, flexural stiffness $\beta$, and unit mass $\alpha$ are positive constants that do not change with respect to $x$ and $y$, such that (2.19) and (2.20) are reduced to

(3.19)$$\begin{gather} \frac{{\partial} \bar{u}_x}{{\partial} t} + \epsilon\bar{u}_x\,\frac{{\partial} \bar{u}_x}{{\partial} x} + \frac{{\partial} \eta}{{\partial} x} + \beta\,\frac{{\partial}^5 \eta}{{\partial} x^5} + \mu^2\alpha\,\frac{{\partial}^3\eta}{{\partial} x\,{\partial} t^2} - \frac{\mu^2 H^2}{3}\,\frac{{\partial}^3 \bar{u}_x}{{\partial} x^2\,{\partial} t} = 0, \end{gather}$$

(3.20)$$\begin{gather}\epsilon\bar{u}_x\,\frac{{\partial} \eta}{{\partial} x} + \epsilon\eta\,\frac{{\partial} \bar{u}_x}{{\partial} x} + H\,\frac{{\partial} \bar{u}_x}{{\partial} x} + \frac{{\partial} \eta}{{\partial} t} = 0. \end{gather}$$

Introducing a moving coordinate $\xi =x-ct$, we have ${ {\partial }}/{ {\partial } x}={ {\partial }}/{ {\partial } \xi }$ and ${ {\partial }}/{ {\partial } t}=-c({ {\partial }}/{ {\partial } \xi })$, such that (3.19) and (3.20) are rewritten as

(3.21)$$\begin{gather} -c\bar{u}_x' + \epsilon\bar{u}_x\bar{u}_x' + \eta' + \beta\eta^{(5)} + \mu^2c^2\alpha\eta''' + \frac{\mu^2 c H^2}{3}\,\bar{u}_x''' = 0, \end{gather}$$

(3.22)$$\begin{gather}\epsilon\bar{u}_x\eta' + \epsilon\eta\bar{u}_x' + H\bar{u}_x' - c\eta' = 0, \end{gather}$$

in which the superscripts indicate the order of derivative with respect to $\xi$.

Equation (3.22) is linearized by assuming $O(\epsilon )\ll O(1)$, which gives $H\bar {u}_x'-c\eta '=0$. Because solitary solutions approach zero for $\xi \rightarrow \pm \infty$, integral constants must be zero and we have $\bar {u}_x=({c}/{H})\eta$, $\bar {u}_x'=({c}/{H})\eta '$ and $\bar {u}_x'''=({c}/{H})\eta '''$. Then (3.21) can be rewritten as a fifth-order nonlinear ordinary differential equation (ODE) with respect to $\eta$:

(3.23)\begin{equation} \left(1-\frac{c^2}{H}\right)\eta' + \epsilon\,\frac{c^2}{H^2}\,\eta\eta' + \mu^2 c^2\left(\alpha+\frac{H}{3}\right)\eta''' + \beta\eta^{(5)} = 0. \end{equation}

It is noted that (3.23) is also equivalent to a generalized fifth-order KdV equation:

(3.24)\begin{equation} \left(\frac{c}{H}-\frac{1}{c}\right)\frac{{\partial}\eta}{{\partial} t} + \epsilon\,\frac{c^2}{H^2}\,\eta\,\frac{{\partial}\eta}{{\partial} x} + \mu^2 c^2 \left(\alpha+\frac{H}{3}\right)\frac{{\partial}^3\eta}{{\partial} x^3} + \beta\frac{{\partial}^5\eta}{{\partial} x^5} = 0. \end{equation}

Substituting a general harmonic solution $\eta =A\exp ({{\rm i}(kx+\sigma t)})$ into the linear part of (3.24), we can derive the approximate dispersion relation of the system: $({(c^2-H)}/{Hc})\sigma =\mu ^2 c^2(\alpha +{H}/{3})k^3+\beta k^5$. The system bifurcates when ${ {\partial } c}/{ {\partial } k}={ {\partial }^2 \sigma }/{ {\partial } k^2}=0$, such that a solitary solution exists at $k\rightarrow 0$ as $c^2\rightarrow H$. Substituting $c^2=H$ into (3.23) and integrating once yields

(3.25)\begin{equation} \eta^{(4)} =-\left(\frac{\mu^2\alpha H}{\beta}+\frac{\mu^2H^2}{3\beta}\right)\eta'' - \frac{\epsilon}{2\beta H}\,\eta^2. \end{equation}

An equivalent dynamical system of (3.25) is derived as

(3.26)\begin{equation} \begin{Bmatrix} \eta \\ \eta' \\ \eta'' \\ \eta''' \end{Bmatrix}'= \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -\left(\dfrac{\mu^2\alpha H}{\beta}+\dfrac{\mu^2 H^2}{3\beta}\right) & 0 \end{bmatrix} \begin{Bmatrix} \eta \\ \eta' \\ \eta'' \\ \eta''' \end{Bmatrix}+ \begin{Bmatrix} 0 \\ 0 \\ 0 \\ -\dfrac{\epsilon\eta^2}{2\beta H} \end{Bmatrix}. \end{equation}

The eigenvalues of the matrix in (3.26) are 0, 0 and $\pm ik_o$, given that $\alpha$, $\beta$, $\epsilon$ and $H$ are positive real numbers. Because the dynamical system has a pair of imaginary eigenvalues, a generalized solitary solution appears, which consists of a ‘${\rm sech}^2$’ core soliton that decays at infinity, and a small-amplitude oscillatory tail with wavenumber $k_o=\sqrt {{\mu ^2\alpha H}/{\beta }+{\mu ^2 H^2}/{3\beta }}$ (Il'ichev Reference Il'ichev1997, Reference Il'ichev2000). The generalized solitary solution for FG waves in shallow water is derived as

(3.27)\begin{equation} \eta_{sol} = \eta_c+\eta_o = 3H\,\frac{c^2-H}{\epsilon c^2}\,{\rm sech}^2 \left[\frac{\tau}{2}\,(x-ct)\right] + A_o\exp\left({-\frac{\rm \pi}{\tau}\,k_o}\right)\exp({{\rm i}k_o(x-ct)}), \end{equation}

in which $\tau =\sqrt {{3(c^2-H)}/{\mu ^2 c^2 H(3\alpha +H)}}$, and $A_o\exp ({-({{\rm \pi} }/{\tau })k_o})$ is the amplitude of the oscillatory tail that is exponentially small with respect to $\sqrt {\beta }$. We present the derivation of the solitary solution $\eta _{sol}$ in Appendix B.

The nonlinear solver developed in § 3.3 is employed to compute the solitary wave given by (3.27), in which the physical parameters from table 1 and water depth 10 m are used. The core soliton $\eta _c$ is used as the initial condition for numerical computations to evaluate the evolution of the solitary solution. The propagating speed $c$ plays an important role in the solitary solution. The solitary solution exists only if $c^2\rightarrow H$. The computation results show that the initial soliton hardly changes when $c$ is close to $\sqrt {H}$ (figure 7a), but breaks down into two major peaks when $c$ gets larger (figure 7b): one peak moves at a lower speed 1.08$\sqrt {H}$ and maintains its shape, while the other propagates at a higher speed 1.32$\sqrt {H}$ and contentiously decays. The amplitude of the core soliton is also dependent on $c$, therefore only small-amplitude solitary waves are permitted.

Figure 7. Solitary waves on sea-ice of constant thickness $d^*=0.1$ m with different wave celerity: (a) $c=1.001\sqrt {H}$, (b) $c=1.1\sqrt {H}$.

Elastic effects are introduced by the fifth-order term in the equivalent KdV equation that generates the oscillatory tail. As demonstrated in Benilov, Grimshaw & Kuznetsova (Reference Benilov, Grimshaw and Kuznetsova1993) and Grimshaw & Joshi (Reference Grimshaw and Joshi1995), the oscillatory tail of a generalized solitary solution is one-sided in a realistic wave system, otherwise a pair of source and sink must be required for energy conservation. The oscillatory tail propagates either in front of or behind the core soliton, depending on the sign of the coefficients in the equivalent KdV equation, and radiates energy from the core. Moreover, a small value of $\beta$ leads to a large wavenumber $k_o$ and results in an exponentially small amplitude of the oscillation with order $\exp [({{-{\rm \pi} c\mu ^2 H(3\alpha +H)}/{3})\sqrt {{1}/{\beta (c^2-H)}}}]$. To evaluate the effects of flexural stiffness $\beta$, solitary waves are simulated on sea-ice with different thicknesses. The numerical results conform to the theoretical predictions. Figure 8 shows that the oscillatory tail appears only in front of the core soliton and has larger wavelength and amplitude for a higher flexural stiffness, resulting in a quicker decay of the core soliton. The simulation of a solitary wave over sea-ice with varying thickness shows similar results (figure 9).

Figure 8. Solitary waves on sea-ices with different thickness (oscillatory tails at $t^*=900$ s are shown in the insets), for $c=1.001\sqrt {H}$ and (a) $d^*=0.6$ m, (b) $d^*=1.2$ m, (c) $d^*=1.8$ m.

Figure 9. Solitary wave propagating over sea-ice with varying thickness (0.1 m to 3 m).

To demonstrate the shoaling effect, a solitary wave propagating in uniform sea-ice ($d^*=0.1$ m) over varying bathymetry is simulated. Figure 10 shows the evolution of the solitary wave when water depth reduces from 10 to 2 m in the range from $x^*=3000$ to $x^*=6000$ m. The speed of the solitary wave reduces, but its amplitude and steepness increase with decreasing water depth. It is expected that the solitary wave would break if it propagates to even shallower water.

Figure 10. Solitary wave on uniform sea-ice over varying water depth ($c=1.001\sqrt {H}$, $d^*=0.1$ m, and $H^*$ changes from $-$10 to $-$2 m).

It should be noted that the solitary solution given by (3.27) is derived assuming that the equivalent KdV equation is singularly perturbed on the fifth-order term, and the KdV equation itself is a reduced model of the proposed Boussinesq equations. Therefore, (3.27) is an approximate solitary solution for FG waves in shallow water with limitations. Exact quantitative analysis of solitary solutions for FG waves in shallow water requires further investigation on the Boussinesq-type equations.

4.1. Dual-mode analytical solutions

The first-order linear governing equations for a submerged plate in shallow water are derived for constant water depth and uniform plate stiffness from (2.28) and (2.29):

(4.1)$$\begin{gather} \frac{{\partial}^2 \hat{\eta}_1}{{\partial} t^2} - \frac{{\partial}^2 \zeta_1}{{\partial} t^2} - H_I\, \frac{{\partial}^2 \hat{\eta}_1}{{\partial} x^2} = 0, \end{gather}$$

(4.2)$$\begin{gather}\frac{{\partial}^2 \zeta_1}{{\partial} t^2} - H_{II}\,\frac{{\partial}^2 \hat{\eta}_1}{{\partial} x^2} - H_{II}\beta\, \frac{{\partial}^6 \zeta_1}{{\partial} x^6} = 0. \end{gather}$$

We solve these equations analytically using separation of variables and periodic conditions by assuming that surface gravity wave $\hat {\eta }_1$ and SH wave $\zeta _1$ have the same periodic characteristics in both space and time, such that $\zeta _1(x,t)=C_3\,\hat {\eta }_1(x,t)=C_3\,\hat {A}(x)\,\hat {B}(t)$. Then (4.1) and (4.2) become

(4.3)$$\begin{gather} \frac{1}{\hat{B}(t)}\,\frac{{\rm d}^2 \hat{B}(t)}{{\rm d} t^2} -\frac{C_3}{\hat{B}(t)}\,\frac{{\rm d}^2 \hat{B}(t)}{{\rm d} t^2} = \frac{H_I}{\hat{A}(x)}\,\frac{{\rm d}^2 \hat{A}(x)}{{{\rm d}x}^2} = C_1, \end{gather}$$

(4.4)$$\begin{gather}\frac{1}{\hat{B}(t)}\,\frac{{\rm d}^2 \hat{B}(t)}{{\rm d} t^2} = \frac{H_{II}}{C_3\hat{A}(x)}\,\frac{{\rm d}^2 \hat{A}(x)}{{{\rm d}x}^2} + \frac{H_{II}\beta}{\hat{A}(x)}\,\frac{{\rm d}^6 \hat{A}(x)}{{{\rm d}x}^6} = C_2, \end{gather}$$

where $C_1$, $C_2$ and $C_3$ are real constants.

Substituting (4.4) into (4.3), we derive a sixth-order homogeneous ODE with respect to $\hat {A}(x)$, whose characteristic equation can be reduced to a cubic equation. Here, $\hat {A}(x)$ is solved by computing the cubic characteristic equation, and the resulting solution for $\hat {A}(x)$ is substituted back into (4.4) to solve $\hat {B}(t)$. Periodic solutions exist in only two scenarios: ${H_{II}}/{(H_I+H_{II})}< C_3<1$ and $C_3<0$. The general periodic solution of $\hat {\eta }_1$ is derived as

(4.5)\begin{equation} \hat{\eta}_1(x,t)=C_4\exp({{\rm i}(\hat{k}x \pm \hat{\sigma} t)}), \end{equation}

in which $\hat {\sigma }$ and $\hat {k}$ are the non-dimensional angular frequency and wavenumber, respectively, and $C_4$ is a real constant determined by boundary and initial conditions.

The SH wave along the submerged plate shares the same frequency and wavenumber as the coupled surface wave as $\zeta _1(x,t)=C_3\,\hat {\eta }_1(x,t)$. Because $C_3$ can be either positive or negative, the coupled hydroelastic wave system has dual-mode analytical periodic solutions: when $C_3$ is positive and always less than $1$, the SH wave has the same phase but a smaller amplitude compared to the coupled gravity wave on the free surface; when $C_3$ is negative, the SH wave is 180$^{\circ }$ out of phase with the surface gravity wave.

In the derivation of analytical solutions, we have

(4.6)\begin{equation} C_3=1-\frac{\hat{k}^2}{\hat{\sigma}^2}\,H_I, \end{equation}

and the dispersion relation in an implicit form is

(4.7)\begin{equation} \hat{k}^2=\sqrt{\frac{\dfrac{\hat{\sigma}^2}{\hat{k}^2}-H_I-H_{II}}{H_{II} \beta-\dfrac{\hat{k}^2}{\hat{\sigma}^2}\,H_IH_{II}\beta}}. \end{equation}

Equation (4.7) can be further reformulated into a quartic equation with respect to $\lambda =\hat {k}^2$:

(4.8)\begin{equation} \frac{H_I H_{II} \beta}{\hat{\sigma}^2}\,\lambda^4 - H_{II}\beta\lambda^3 - (H_I+H_{II})\lambda + \hat{\sigma}^2=0. \end{equation}

Solving the quartic equation with Ferrari's method, (4.8) has two positive real roots and a pair of complex conjugate roots, given that $H_I$, $H_{II}$, $\beta$ and $\hat {\sigma }$ are positive real constants. The wavenumbers of the periodic solutions are the square roots of the two positive real roots of $\lambda$, corresponding to ${H_{II}}/{(H_I+H_{II})}< C_3<1$ and $C_3<0$, respectively. Considering (4.6), the smaller root of $\lambda$ gives the wavenumber of the in-phase mode with larger wavelength and celerity $({L}/{T})^2>H_I+H_{II}$, and the greater root of $\lambda$ yields the wavenumber of the out-of-phase mode with smaller wavelength and celerity $({L}/{T})^2< H_I$. Hereafter, we employ superscripts to denote the wavenumber and $C_3$: $\hat {k}^+$ and $C_3^+$ are associated with the in-phase mode, and $\hat {k}^-$ and $C_3^-$ are associated with the out-of-phase mode.

To demonstrate the dual-mode solutions and the effects of submerged depth, we compute the gravity–hydroelastic coupled waves on a submerged polyvinyl chloride (PVC) plastic plate, using the analytical solutions. Table 3 presents the physical parameters used for the computations. Computing and substituting the corresponding dimensionless parameters ($H_I$, $H_{II}$, $\beta$ and $\hat {\sigma }$) into (4.8) to solve for $\lambda$ and $\hat {k}$, the dual-mode solutions of $\hat {\eta }_1$ and $\zeta _1$ are obtained using (4.5) and (4.6). Figures 11 and 12 show the in-phase and out-of-phase coupled waves computed for various submerged depths, respectively.

Table 3. Physical parameters for submerged PVC plate computation.

Figure 11. Analytical solutions (in-phase mode) with a PVC plate submerged at different depths (wave amplitudes are exaggerated by a factor of 10). Submerged depths are (a) $1$ m, (b) $2$ m, (c) $3$ m.

Figure 12. Analytical solutions (out-of-phase mode) with a PVC plate submerged at different depths (wave amplitudes are exaggerated by a factor of 10). Submerged depths are (a) $1$ m, (b) $2$ m, (c) $3$ m.

For the in-phase mode (figure 11), the amplitude of the SH wave decreases with increasing submerged depth, while wavelength and celerity of the coupled waves decrease with increasing submerged depth. For the out-of-phase mode (figure 12), the wavelength, wave celerity and amplitude of the SH wave all increase with increasing submerged depth. Amplitude of the SH wave can be even greater than that of the surface wave in the out-of-phase mode when $C_3<-1$ (see figure 12c).

4.2. First-order linear dispersion relation

To further investigate the characteristics of the dual-mode solutions, we rewrite (4.7) and (4.8) for the first-order linear dispersion relation:

(4.9)\begin{equation} \frac{H_I H_{II} \beta}{\hat{\sigma}^2}\,\hat{k}^8 - H_{II}\beta \hat{k}^6 - (H_I+H_{II})\hat{k}^2 + \hat{\sigma}^2=0. \end{equation}

Equation (4.9) reduces to (3.2), the first-order dispersion relation of FG waves, when the submerged depth $H_I$ goes to zero. Equation (4.9) also reduces to the first-order dispersion relation of shallow-water waves when $H_{II}$ goes to zero ($H_I\rightarrow H$) or $\beta$ becomes zero. The linear dispersion relation for SH waves in water of finite depth was derived by Hassan et al. (Reference Hassan, Meylan and Peter2005), which can be non-dimensionalized as

(4.10)\begin{align} & \mu \hat{\sigma}^2 \beta \hat{k}^5 - \beta \hat{k}^6\tanh(\mu \hat{k}H_I) - \mu^3 \alpha \hat{\sigma}^4\hat{k} + \mu^2\alpha\hat{\sigma}^2\hat{k}^2\tanh(\mu \hat{k}H_I) \nonumber\\ &\quad =\mu^2\hat{\sigma}^4\tanh(\mu \hat{k}H_I)-\mu\hat{\sigma}^2\hat{k} - \frac{\mu\hat{\sigma}^2\hat{k}\tanh(\mu \hat{k}H_I)-\mu^2\hat{\sigma}^4}{\tanh(\mu \hat{k}H_{II})}. \end{align}

Equation (4.10) reduces to (4.9) for shallow water, when $\mu$ is small ($\tanh \mu \hat {k}H\rightarrow \mu \hat {k}H$) and its higher-order terms are neglected.

The effects of submerged depth on SH waves are investigated for various wave frequencies and total depths. The dispersion relation (4.9) is solved for $\hat {k}^+$ and $\hat {k}^-$ for different submerged depths $H_I$, and the resulting $C_3^+$ and $C_3^-$ are computed as well. Figures 13 and 14 show that the wavelength and $C_3^+$ decrease with increasing submerged depth for the in-phase mode. For the out-of-phase mode, on the other hand, its wavelength generally increases with increasing submerged depth. Such results can be explained by combining (2.28) and (2.29), which yields

(4.11)\begin{equation} \frac{{\partial}^2 \hat{\eta}_1}{{\partial} t^2} - \boldsymbol{\nabla}\boldsymbol{\cdot}(H\,\boldsymbol{\nabla}\hat{\eta}_1) - \boldsymbol{\nabla}\boldsymbol{\cdot}\{H_{II}\,\boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\zeta_1)]\} =0. \end{equation}

The first two terms in (4.11) are similar to those in the linear FG wave equation (2.17), representing the gravitational effects of surface waves, and the third term represents the elastic effects of the submerged plate. Comparing (4.11) with (2.17), it is noted that the effective depth of the elastic term becomes the depth below the submerged plate ($H_{II}$), such that the wavelength decreases with decreasing $H_{II}$ (increasing $H_I$) for the in-phase mode. Substituting the out-of-phase mode $\zeta _1=-|C_3|\,\hat {\eta }_1$ into the third term of (4.11), the sign of the elastic term flips comparing to the in-phase mode, resulting in the opposite influence of submerged depth on elastic effects.

Figure 13. Dispersion relations of SH waves in shallow water with various wave periods ($L_0=50$ m, $H_0 =4\ {\rm m}$, $\beta =3.155\times 10^{-4}$, PVC thickness 0.4 m).

Figure 14. Dispersion relations of SH waves in shallow water with various water depths ($L_0=50$ m, $T^{*}=5\ {\rm s}$, $\beta =3.155\times 10^{-4}$, PVC thickness 0.4 m).

The effects of flexural stiffness on the first-order dispersion relation are investigated similarly. Figure 15 shows that increasing plate stiffness results in increase of wavelength for both in-phase and out-of-phase modes. When $\beta \rightarrow 0$, the elastic effects fade away and the in-phase mode converges to the normal shallow-water wave, such that $\hat {k}^+$ becomes independent of submerged depth, and $C_3^+$ presents motion magnitude of a submerged water particle relative to the surface. When $\beta \rightarrow \infty$, the submerged plate plays a role similar to the rigid bottom boundary, and the out-of-phase mode converges to the normal shallow-water wave over the rigid plate, such that $\hat {k}^-$ approaches the wavenumber of the shallow-water wave, and $C_3^-$ reaches zero.

Figure 15. Dispersion relations of SH waves in shallow water with various flexural stiffnesses ($L_0=50\ {\rm m}$, $H_0=4$ m, $T^*=5$ s).

4.3. Dual-mode numerical solutions

We develop a numerical solver to compute the first-order gravity–hydroelastic coupled wave system, which allows for demonstration of the dual-mode solutions with varying flexural stiffness and varying submerged depth.

The equations preceding (2.28) and (2.29), other than themselves, are solved to reduce the order of time derivatives:

(4.12)$$\begin{gather} \frac{{\partial} \hat{\eta}_1}{{\partial} t} - \frac{{\partial} \zeta_1}{{\partial} t} =- \boldsymbol{\nabla}\boldsymbol{\cdot}(H_I\bar{\boldsymbol{u}}_{I1}), \end{gather}$$

(4.13)$$\begin{gather}\frac{{\partial} \zeta_1}{{\partial} t} =- \boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1}), \end{gather}$$

(4.14)$$\begin{gather}\frac{{\partial} \bar{\boldsymbol{u}}_{I1}}{{\partial} t} =- \boldsymbol{\nabla}\hat{\eta}_1, \end{gather}$$

(4.15)$$\begin{gather}\frac{{\partial} \bar{\boldsymbol{u}}_{II1}}{{\partial} t} =- \boldsymbol{\nabla}\hat{\eta}_1 - \boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\zeta_1)]\}, \end{gather}$$

in which $\bar {\boldsymbol {u}}_{I1}$ and $\bar {\boldsymbol {u}}_{II1}$ are the first-order horizontal velocities above and below the submerged plate, respectively. The use of provisional values allows conditionally stable computation using the central difference schemes for spatial derivatives. Therefore, provisional steps are used in an explicit time-marching scheme to solve (4.12)–(4.15):

(4.16)$$\begin{gather} \zeta_{1,s+1^*} = \zeta_{1,s} - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1,s}), \end{gather}$$

(4.17)$$\begin{gather}\hat{\eta}_{1,s+1^*} = \hat{\eta}_{1,s} - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1,s}) - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{I}\bar{\boldsymbol{u}}_{I1,s}), \end{gather}$$

(4.18)$$\begin{gather}\bar{\boldsymbol{u}}_{I1,s+1} = \bar{\boldsymbol{u}}_{I1,s} - \Delta t\,\boldsymbol{\nabla}\left(\frac{\hat{\eta}_{1,s}+\hat{\eta}_{1,s+1^*}}{2}\right)\!, \end{gather}$$

(4.19)$$\begin{gather}\bar{\boldsymbol{u}}_{II1,s+1} = \bar{\boldsymbol{u}}_{II1,s} - \Delta t\,\boldsymbol{\nabla}\left(\frac{\hat{\eta}_{1,s}+\hat{\eta}_{1,s+1^*}}{2}\right) - \Delta t\,\boldsymbol{\nabla}\left\{\nabla^2\left[\beta\,\nabla^2\left(\frac{\zeta_{1,s}+\zeta_{1,s+1^*}}{2}\right)\right]\right\}\!, \end{gather}$$

(4.20)$$\begin{gather}\zeta_{1,s+1} = \zeta_{1,s} - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1,s+1}), \end{gather}$$

(4.21)$$\begin{gather}\hat{\eta}_{1,s+1} = \hat{\eta}_{1,s} - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1,s+1}) - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{I}\bar{\boldsymbol{u}}_{I1,s+1}), \end{gather}$$

in which $\hat {\eta }_{1,s+1^*}$ and $\zeta _{1,s+1^*}$ are the provisional values of $\hat {\eta }_{1,s+1}$ and $\zeta _{1,s+1}$, respectively. A linear numerical solver is then developed based on (4.16)–(4.21).

The linear gravity–hydroelastic coupled wave system is computed numerically to demonstrate the influence of varying flexural stiffness. A harmonically oscillating boundary excitation is applied to the left end on the free surface to induce a regular incident wave. Figure 16 presents the responses of the free surface and the PVC plates with different stiffnesses, in which the parameters in table 3 are employed for computations, except the thickness. Figure 17 shows the coupled wave responses with plate thickness varying from $d^*=0.1$ to $d^*=0.6$ m. It is clear that both the in-phase and out-of-phase modes co-exist. The surface wave does not change noticeably due to the varying stiffness, while the SH wave is significantly affected by the plate stiffness: the out-of-phase mode is more prominent in the plate with lower stiffness, which conforms to the observation in figure 15(d).

Figure 16. Linear gravity–hydroelastic coupled waves with different flexural stiffnesses ($H_0=4$ m, $H_I=0.5$, wave amplitudes are exaggerated by a factor of 10). (a) $\beta = 1.204 \times 10^{-5}\ (d^{\ast} = 0.1\ {\rm m})$; (b) $\beta = 3.250 \times 10^{-4} (d^{\ast} = 0.3\ {\rm m})$; (c) $\beta = 2.600 \times 10^{-3}\ (d^{\ast} = 0.6\ {\rm m})$.

Figure 17. Linear gravity–hydroelastic coupled waves with varying flexural stiffnesses ($H_0=4$ m, $H_I=0.5$, wave amplitudes are exaggerated by a factor of 10). (a) $\beta{:}\ 1.204 \times 10^{-5}\unicode{x2013}2.600 \times 10^{-3}$; (b) $\beta{:}\ 2.600 \times 10^{-3}\unicode{x2013}1.204 \times 10^{-5}$.

The influence of varying submerged depth on the coupled wave system is more complicated than that of varying stiffness. Figure 18 shows the responses of the free surface and the PVC plate submerged at different depths subject to the harmonic boundary excitation. Both the surface elevation and plate motion are sensitive to submerged depth: the dominant mode of the coupled wave system changes from the in-phase mode to the out-of-phase mode with increasing submerged depth. Figure 19 presents the coupled wave responses with submerged depth varying from 0.5 to 3.5 m in total water depth of 4 m. When submerged depth gradually increases from $H_I=0.125$ to $H_I=0.875$ (figure 19a), the in-phase mode is dominant: the amplitude of the SH wave decreases dramatically with increasing submerged depth, while the amplitude of the surface wave sightly increases, which conforms to figures 13(b) and 14(b). When submerged depth gradually decreases from $H_I=0.875$ to $H_I=0.125$ (figure 19b), the out-of-phase mode component in the surface wave significantly increases, while the in-phase mode component in the SH wave increases.

Figure 18. Linear gravity–hydroelastic coupled waves with various submerged depths ($d^*=0.2$ m, $\beta =9.629\times 10^{-5}$, wave amplitudes are exaggerated by a factor of 10). Submerged depths are (a) $1$ m, (b) $2$ m, (c) $3$ m.

Figure 19. Linear gravity-hydroelastic-coupled waves with varying submerged depth ($d^*=0.4$ m, $\beta =3.155\times 10^{-4}$, wave amplitudes are exaggerated by a factor of 10). (a) $H_I : 0.125\unicode{x2013}0.875$; (b) $H_I : 0.875\unicode{x2013}0.125$.

Energy flows into the coupled wave system from the boundary excitation and is distributed between the in-phase and out-of-phase modes. For each mode, the modal energy is further distributed between the free surface and the submerged plate based on the values of $C_3^+$ and $C_3^-$. The coupled wave responses are the combined results of these energy distributions. The result in figure 19(b) is appealing because we notice that the surface wave energy density remarkably increases due to the plate with decreasing submerged depth. The energy input from the oscillating boundary is allocated to the kinetic energies of the upper and lower flows, the gravitational potential at the free surface, and the elastic and kinetic energies of the submerged plate. The kinetic energy associated with the plate is neglected in the linear model, while the elastic potential is included. The free surface wave is associated with the gravitation potential and the kinetic energy of the upper flow, and similarly the SH wave is an indicator of the elastic potential of the plate and the kinetic energy of the lower flow. Figure 20 shows the horizontal velocities of the upper and lower flows. The changes in flow velocities follow the coupled waves in figure 19(b): the amplitude of the upper flow velocity increases, and its dominant mode becomes the out-of-phase mode as the submerged depth decreases. The average energy per unit surface area over a wave period can be computed based on the formula

(4.22)\begin{align} P(x^*,t^*) &= \frac{1}{2T^*}\,\rho \int_{t^*}^{t^*+T^*} \int_{-H_I^*+\zeta_1^*}^{\hat{\eta}_1^*}\bar{u}_{I1}^{*2} \,{\rm d} z^*\,{\rm d} t^* + \frac{1}{T^*}\,\rho \int_{t^*}^{t^*+T^*}\int_{0}^{\hat{\eta}_1^*} gz^*\,{\rm d} z^*\,{\rm d} t^* \nonumber\\ &\approx \frac{1}{2T^*}\,\rho H_I^* \int_{t^*}^{t^*+T^*} \bar{u}_{I1}^{*2} \,{\rm d} t^* + \frac{1}{2T^*}\,\rho g \int_{t^*}^{t^*+T^*} \hat{\eta}_1^{*2} \,{\rm d} t^*. \end{align}

Both the in-phase and out-of-phase modes share the same period as the boundary excitation ($T^*=10$ s in these computations). Substituting the time histories of $\bar {u}_{I1}^*$ and $\hat {\eta }_1^*$ into (4.22), we compute that the surface wave energy density at $x^*=600$ m is nearly nine times greater than that at $x^*=100$ m (7.86 J m$^{-2}$ versus 0.79 J m$^{-2}$). Therefore, the plate with decreasing submerged depth effectively extracts wave energy from the deeper water and transfers it to the surface layer. Such effect is enabled by energy redistributions between the in-phase and out-of-phase modes, and between the surface and SH waves. We have investigated the energy distribution between the free surface and SH waves based on the values of $C_3^+$ and $C_3^-$ for various submerged depth and flexural stiffness. However, the energy distribution between the in-phase and out-of-phase modes is not trivial. That depends on total water depth, submerged depth, properties of the submerged plate, and the means of energy input, such as frequency and location of the excitation. Analysis on the dual-mode solutions of SH waves will be within the scope of our future study.

Figure 20. Averaged horizontal velocities of the upper and lower flows with varying submerged depths (3.5–0.5 m).