Analytical Solution for Wave Scattering by a Surface Obstacle above a Muddy Seabed

Abstract

We present an analytical solution for the scattering of linear progressive waves by a surface rectangular obstacle above a muddy seabed. The bottom cohesive mud is assumed to act as a Newtonian fluid, and the thickness of the mud layer is considered to be comparable to the Stokes boundary layer thickness. Our analytical results based on the matched eigenfunction expansions incorporate the combined effects of obstacles and a fluid mud bottom. By reducing the mud layer thickness or the dimensions of the obstacle to zero, the present study recovers the classical solution for wave scattering by a surface obstacle above a solid bed or wave propagation over a layer of fluid mud. Our analytical predictions of wave amplitudes and wave forces acting on the bottom of the obstacle agree satisfactorily with the available numerical results. The most prominent effect of a muddy seabed is a strong damping of wave amplitude. Parameter study reveals that the obstacle submerged depth, mud layer thickness, and wave frequency can have significant impacts on the attenuation of wave amplitude due to the presence of a muddy seabed.

1. Introduction

The scattering of surface water waves by a fixed structure of finite dimensions remains an intruding and practical research topic in the engineered coastal and offshore systems. Many classical analytical solutions based on the linear wave theory have been reported for waves interacting with a floating or submerged obstacle [1,2,3,4,5]. For instance, Mei and Black [4] employed the eigenfunction expansions [1] to study the effects of a fixed and rigid rectangular prism on two-dimensional wave propagation on a constant water depth. Prominent scattering features such as reflection and transmission coefficients were reported. The method of eigenfunction expansions is one of the well-received mathematical techniques for wave-structure problems. It has been applied to many related problems, such as diffraction of obliquely incident waves by a rectangular trench [6], wave scattering by a porous structure [7], and even the wave problem involving a moving obstacle with small displacements [8]. More recent analytical studies on wave-obstacle problems include the use of mild-slope assumption to investigate the effects of bottom undulation on wave scattering [9], the application of dual poroviscoelastic wave barriers as a protection measure for floating solar farms [10], and the scattering of flexural-gravity waves by a floating ice sheet [11], among others.

In many coastal areas where the obstacles of interest can be relevant for various engineering applications, bottom cohesive sediments may interact with the wave field. Consequently, the wave–seabed interaction not only changes wave characteristics such as wave celerity, wavelength, and wave height but also causes significant motions of cohesive sediments [12]. The effects of a muddy seabed on wave propagation have been frequently observed in the field [13,14]. Past studies of wave–seabed interactions often employed a two-layer system consisting of a upper inviscid water layer and a bottom muddy seabed [12,15,16,17,18]. Suggested by the rheological testing of field samples taken from various sites, several representative theoretical rheological models, including Newtonian fluids [12,16], viscoelastic fluids [15,17,18], viscoplastic materials [19], and poroelastic models [20], have been proposed to model the rheological behaviors of the bottom fluid mud. Due to the diverse rheological properties of fluid mud partly as a consequence of its distinct physico-chemical compositions, it is unequivocal that no single model is capable of describing the entire spectrum of rheological behaviors of fluid mud. Recent analytical results for wave–mud interactions include a depth-integrated model for weakly nonlinear long waves over a thin layer of viscoelastic mud on a mild-slope beach [21], a Boussinesq-type model incorporating two distinct soft mud layers to take into account the vertical variation of mud properties [22], and a linear theory for waves and currents over a viscoud fluid bed [23].

Through the above discussion, it is evident that continuous efforts have been made to study both wave–obstacle and wave–mud interactions. However, the combined wave–obstacle–mud problem seems to receive less attention. Therefore, in the present study, we shall investigate the scattering of water waves by a surface obstacle above a muddy seabed. A common two-layer system is adopted for the wave–obstacle–seabed problem. The bottom fluid mud is idealized as a Newtonian fluid. We consider only waves of small amplitudes in a two-dimensional setting. The method of eigenfunction expansions is employed as the analytical tool to incorporate the effects of a surface obstacle in the shape of a rectangle. Due to the presence of a muddy seabed, the usual orthogonality of eigenfunctions can no longer be applied to deduce unknown coefficients in the expansions. In this regard, the idea of orthogonal mode-coupling [11,24,25] is adopted.

In Section 2, we present the formulations and the solution forms for wave scattering by a surface obstacle above a muddy seabed. Assumptions and simplifications are also discussed. Section 3 is devoted to the solution methods for the unknown coefficients and the dispersion relations. Model predictions are examined using a set of available numerical results [26], and a reasonable agreement is observed as presented in Section 4. With the help of the new solution, the effects of key parameters, such as wave frequency, obstacle submerged depth, and mud layer thickness, on wave amplitudes and wave forces are discussed and presented also in Section 4. Key findings of the present study and the limitations of the present model are summarized and discussed in Section 5. Finally, we conclude with the potential directions for future studies.

2. Solution Forms

2.1. Assumptions and Simplifications

We consider a two-dimensional two-layer system commonly used in the theoretical studies of wave–mud interactions [12], as sketched in Figure 1. The upper water layer follows the typical assumptions of a perfect fluid and linear progressive waves [4]. The bottom muddy seabed is idealized by a Newtonian fluid. Although the simple viscous model has been shown to efficiently capture the prominent features of mud-induced wave energy dissipation [14], we reiterate that there are still many occasions in which the Newtonian model fails to accurately describe the field reality due to the complex rheological properties of fluid mud, as discussed previously.

Assuming laminar flows and the thickness of the mud layer to be much smaller than the water depth but comparable to the Stokes boundary layer thickness of mud, interfacial mixing can be ignored, and a sharp interface is persistent between the two layers [18]. A solid obstacle in the shape of a rectangle is fixed in the free surface. The surface obstacle is sufficiently high so that no wave overtopping is possible.

2.2. Muddy Seabed

Inside the muddy seabed, the motion of the fluid mud is governed by the linearized Navier–Stokes equations [15],

$$\frac{\partial u_m}{\partial x}+\frac{\partial w_m}{\partial z}=0,$$

(1)

$$\frac{\partial u_m}{\partial t}=-\frac{1}{{\rho }_m}\frac{\partial p_m}{\partial x}+{\nu }_m\left(\frac{{\partial }^2{u}_m}{\partial x^2}+\frac{{\partial }^2{u}_m}{\partial z^2}\right),$$

(2)

$$\frac{\partial w_m}{\partial t}=-\frac{1}{{\rho }_m}\frac{\partial p_m}{\partial z}+{\nu }_m\left(\frac{{\partial }^2{w}_m}{\partial x^2}+\frac{{\partial }^2{w}_m}{\partial z^2}\right),$$

(3)

where t is time, $(x,z)$ are defined in Figure 1, $(u_m,w_m)$ are velocity components, $p_m$ is pressure, ${\rho }_m$ is mud density, and ${\nu }_m$ is the viscosity of mud. We reiterate that nonlinear convection terms have been neglected, since we focus only on a linear problem. At the very bottom, the no-slip condition requires

$$u_m=0,z=-(h+d),$$

(4)

and

$$w_m=0,z=-(h+d),$$

(5)

where h is the constant depth in the incident region and d is the thickness of the mud layer. At the water–mud interface, the dynamic conditions require the continuity of tangential and normal stress components, respectively, as [18]

$${\rho }_m{\nu }_m\frac{\partial u_m}{\partial z}=0,z=-h,$$

(6)

and

$$-p_m+2{\rho }_m{\nu }_m\frac{\partial u_m}{\partial z}=-p_w,z=-h,$$

(7)

where $p_w$ is the pressure of water above the muddy seabed. In addition, the kinematic condition states that

$$w_w=w_m=\frac{\partial {\zeta }_m}{\partial t},z=-h,$$

(8)

where $w_w$ is the vertical component of water particle velocity and ${\zeta }_m(x,t)$ denotes the interfacial displacement.

We assume that d is much smaller than the water depths but comparable to the Stokes boundary layer thickness defined as

$${\delta }_m=\sqrt{\frac{2{\nu }_m}{\omega }},$$

(9)

where $\omega$ is the wave frequency. It follows that the horizontal momentum Equation (2), can be approximated by the typical boundary layer equation as [17]

$$\frac{\partial u_m}{\partial t}=-\frac{1}{{\rho }_m}\frac{\partial p_m}{\partial x}+{\nu }_m\frac{{\partial }^2{u}_m}{\partial z^2}$$

(10)

and the vertical momentum Equation (3), suggests a vertically uniform pressure inside the mud layer [17],

$$p_m=p_m(x,t).$$

(11)

Since in the present study $\mathcal{O}\left(d\right)=\mathcal{O}\left({\delta }_m\right)$, the dynamic condition in the normal direction, (7), can be simplified as the continuity of pressure across the interface [17],

$$p_m=p_w,z=-h.$$

(12)

Hence, using (12) in (10), the horizontal momentum Equation (10), becomes

$$\frac{\partial u_m}{\partial t}=-\frac{1}{{\rho }_m}\frac{\partial p_b}{\partial x}+{\nu }_m\frac{{\partial }^2{u}_m}{\partial z^2},$$

(13)

where

$$p_b=p_w(x,z=-(h+d),t).$$

(14)

Flow motion in the water layer is governed by the Euler equations. Therefore, the pressure gradient term in (13) can be expressed as

$$\frac{\partial p_b}{\partial x}=-{\rho }_w\frac{\partial u_w}{\partial t},$$

(15)

where $u_w$ is the horizontal water particle velocity and ${\rho }_w$ is the density of water. Substituting (15) into (13), the linearized horizontal momentum equation in the muddy seabed becomes [18]

$$\frac{\partial u_m}{\partial t}=\frac{{\rho }_w}{{\rho }_m}\frac{\partial u_b}{\partial t}+{\nu }_m\frac{{\partial }^2{u}_m}{\partial z^2},$$

(16)

where

$$u_b=u_w(x,z=-(h+d),t).$$

(17)

The two-point boundary value problem for $u_m$, i.e., (16) with (4) and (6) can be readily solved to obtain the solution form for $u_m$ as [18]

$$u_m(x,z,t)=\frac{{\rho }_w}{{\rho }_m}{u}_b\left[-cosh\alpha \left(z+h+d\right)+tanh\left(\alpha d\right)sinh\alpha \left(z+h+d\right)+1\right],$$

(18)

where

$$\alpha =\left(1-\mathrm{i}\right)\sqrt{\frac{\omega }{2{\nu }_m}}.$$

(19)

Consequently, by the conservation of mass, (1), the solution for $w_m$ is

$$\begin{array}{cc}\hfill w_m(x,z,t)=& {\int }_{-(h+d)}^z-\frac{\partial u_m(x,z^{\prime },t)}{\partial x}d{z}^{\prime }\hfill \\ \hfill =& \frac{{\rho }_w}{{\rho }_m}\frac{\partial u_b}{\partial x}\left\{\frac{sinh\alpha \left(z+h+d\right)}{\alpha }\right.\hfill \\ \hfill & \left.-\frac{tanh\alpha d}{\alpha }\left[cosh\alpha \left(z+h+d\right)-1\right]-\left(z+h+d\right)\right\}.\hfill \end{array}$$

(20)

We note that in (18) and (20), the water particle velocity $u_b(x,t)$, defined in (17), is yet to be determined. For later use, the vertical velocity at the water–mud interface is

$$w_b=w_w{|}_{z=-h}=w_m{|}_{z=-h}=\frac{{\rho }_w}{{\rho }_m}\frac{\partial u_b}{\partial x}\left[d-\frac{tanh\alpha d}{\alpha }\right].$$

(21)

2.3. Water Layer

With the usual assumptions for water waves, we can define a velocity potential for the motion of water particles as

$$\nabla \Phi =(u_w,w_w),$$

(22)

which satisfies the Laplace equation

$${\nabla }^2\Phi =\frac{{\partial }^2\Phi }{\partial x^2}+\frac{{\partial }^2\Phi }{\partial z^2}=0,\left\{\begin{array}{c}\mathrm{Region}1:x<-B,-h<z<0,\hfill \\ \mathrm{Region}2:-B<x<B,-H<z<-W,\hfill \\ \mathrm{Region}3:B<x,-h<z<0,\hfill \end{array}\right.$$

(23)

and the following conditions

$$\left\{\begin{array}{ccc}\frac{\partial \Phi }{\partial z}=\frac{\partial {\zeta }_w}{\partial t},\frac{\partial \Phi }{\partial t}=-g{\zeta }_w,z=0\hfill & \hfill & (\mathrm{Region}1\&\mathrm{Region}3:|x|>B)\hfill \\ \frac{\partial \Phi }{\partial z}=0,z=-W\hfill & \hfill & \left(\mathrm{Region}2:\right|x|<B)\hfill \end{array}\right.,$$

(24)

$$\frac{\partial \Phi }{\partial z}=w_m,z=0,$$

(25)

$$\frac{\partial \Phi }{\partial x}=0,\left\{\begin{array}{c}x=\pm B\hfill \\ -W<z<0\hfill \end{array}\right..$$

(26)

In the above, g is the gravitational acceleration, ${\zeta }_w$ is the free-surface displacement, H is the constant depth below the surface obstacle of length $2B$ and submerged depth W. In addition, the pressure filed can be obtained from the unsteady Bernoulli equation as

$$p_w=\left(\mathrm{hydrostatic}\right)+\left(\mathrm{dynamic}\right)=-{\rho }_{w}gz+{\rho }_w\frac{\partial \Phi }{\partial t}.$$

(27)

We consider first the flow region before the surface obstacle, i.e., Region 1: $x<-B$. To satisfy the Laplace Equation (23), by the method of eigenfunction expansions, the solution form for $\Phi$ is formulated as

$$\begin{array}{cc}\hfill {\Phi }^{\left(1\right)}=e^{-\mathrm{i}\omega t}\{& A_1{e}^{\mathrm{i}{k}_1(x+B)}\left[cosh{k}_1(z+h)+{\tilde{\Omega }}_{1}sinh{k}_1(z+h)\right]\hfill \\ \hfill & +\sum _{n=1}^{\infty }{B}_n{e}^{-\mathrm{i}{k}_n(x+B)}\left[cosh{k}_n(z+h)+{\tilde{\Omega }}_{n}sinh{k}_n(z+h)\right]\}\hfill \end{array},$$

(28)

where $(A_1,B_n,{\tilde{\Omega }}_n)$ and $k_n$ are unknowns to be determined and the infinite terms appear due to the presence of the surface obstacle [4]. We remark that the term associated with $A_1$ is related to the incident waves while $B_1$ term represents the reflected waves [4]. For convenience, the superscript (1) is also introduced to denote Region 1. Substituting (28) into the interfacial condition, (25), and using (21), we obtain the coefficient ${\tilde{\Omega }}_n$ to be

$${\tilde{\Omega }}_n=\Omega k_n,\Omega =\frac{{\rho }_w}{{\rho }_m}\left(\frac{tanh\alpha d}{\alpha }-d\right).$$

(29)

It follows from the free-surface conditions, (24), that a dispersion relation can be derived as [18]

$${\omega }^2=g{k}_{n}tanh{k}_{n}h\frac{1-\Omega coth{k}_{n}h}{1-\Omega tanh{k}_{n}h}.$$

(30)

Clearly, $\Omega$ represents the effects of a muddy seabed on wave propagation. The dispersion relation needs to be solved numerically and is discussed shortly in Section 3.1.

Similarly, in Region 3, i.e., $x>B$, the solution form for the velocity potential can be readily obtained as

$${\Phi }^{\left(3\right)}=e^{-\mathrm{i}\omega t}\sum _{n=1}^{\infty }{E}_n{e}^{\mathrm{i}{k}_n(x-B)}\left[cosh{k}_n(z+h)+\Omega k_{n}sinh{k}_n(z+h)\right],$$

(31)

where $E_n$ is to be determined and $\Omega$ and $k_n$ are the same as those given in (29) and (30), respectively. The flow motion is expected to be different in Region 2 $\left(\right|x|<B)$, since the bottom of the surface solid obstacle acts as a rigid lid, as indicated by the no-flux condition in (24). Following the same procedure, we shall obtain the velocity potential for $\left|x\right|<B$ as [4]

$${\Phi }^{\left(2\right)}=e^{-\mathrm{i}\omega t}\left\{C_1+D_1\frac{x}{B}+\sum _{m=2}^{\infty }\left[C_m\frac{cosh{K}_{m}x}{cosh{K}_{m}B}+D_m\frac{sinh{K}_{m}x}{sinh{K}_{m}B}\right]cos{K}_m(z+W)\right\},$$

(32)

where $(C_m,D_m)$ are unknown coefficients and $K_m$ satisfies

$$\Omega K_m=tan{K}_{m}H.$$

(33)

We have obtained the complete solution form for the velocity potential. In the following, we first comment on the solution procedure for the dispersion relations, (30) and (33); then, we discuss the approach to obtain the unknown coefficients, $(A_n,B_n,C_m,D_m,E_n)$.

3. Dispersion Relations and Unknown Coefficients

3.1. Dispersion Relations

The dispersion relation for both Region 1 and Region 3 is given in (30). Clearly, in the absence of a muddy seabed, the mud response function, $\Omega$ defined in (29), reduces to $\Omega =0$ and, consequently, (30) recovers the familiar form for the case of a solid bottom with $k_1$ being the wavenumber of the incident wave, ${\tilde{k}}_0$, and $k_{n>1}$ representing the evanescent mode [4]. If there is no surface obstacle, (30) permits only one relevant solution, a complex wavenumber $\tilde{k}$ whose real part represents wave propagation, and the imaginary part shows the wave amplitude attenuation rate due to bottom mud [17]. Furthermore, if the mud thickness is much smaller than the typical wavelength of the incident wave, i.e., $\epsilon =\mathcal{O}({\tilde{k}}_{0}d)\ll 1$, an asymptotic solution for $\tilde{k}$ can be obtained by the perturbation expansion in $\epsilon$ as [18]

$$\tilde{k}\approx \mathcal{O}\left(1\right)+\mathcal{O}\left(\epsilon \right)={\tilde{k}}_0+{\tilde{k}}_1,\left\{\begin{array}{c}{\omega }^2=g{\tilde{k}}_{0}tanh{\tilde{k}}_{0}h\hfill \\ {\tilde{k}}_1=\frac{\Omega {\tilde{k}}_0}{sinh{\tilde{k}}_{0}hcosh{\tilde{k}}_{0}h+{\tilde{k}}_{0}h}\hfill \end{array}\right..$$

(34)

In the present study, the approximation given in (34) is used as the initial guess for solving numerically the full dispersion relation (30). We remark that an infinite number of $k_n$ shall be obtained. However, only the first several terms are important [4]. The numerical implement for solving (30) is standard. For instance, the common fsolve syntax in either MATLAB or Python is sufficient for the computational task. We do not reiterate here.

As an example, Figure 2 demonstrates the wavenumbers for the case of $h=10$ m, ${\rho }_w/{\rho }_m=0.9$, ${\nu }_m=0.01$ m${}^2$/s, and $d/{\delta }_m=1.5$. As shown in the figure, all $k_n$ are complex numbers. The imaginary part of the propagation mode, see $k_1$ in the top panel of the figure, represents the wave amplitude attenuation rate due to the presence of a muddy seabed [17]. In addition, stronger damping occurs in shallow waters. As for the solutions of $k_{n>1}$ shown in the bottom panel of Figure 2, they represent a rapid exponential decay away from the obstacle, which is typical for wave-structure problems. As can be expected, the real part of $k_1$ is much larger than those of $k_{n>1}$ while the opposite is true for the imaginary parts. We remark that for the evanescent mode, $k_{n>1}$, only the imaginary parts are relevant.

Regrading the dispersion relation for Region 2, since (33) is in the form similar to the typical equation governing the evanescent mode [4], the numerical solution for $K_m$ follows the standard procedure given in the literature [27]. We note that by setting $d=0$, (33) reduces to the case of a solid seabed and $K_m=(m-1)\pi /H$ [4]. Figure 3 shows the results of $K_{m>1}$ under $H/h=0.5$. In the present study, we have complex wavenumbers, and the real part of $K_m$ is approximately $\Re \left(K_m\right)=(m-1)\pi /H$.

3.2. Unknown Coefficients

If the free-surface displacement of the incident wave is formulated as

$${\zeta }_w=a_0{e}^{\mathrm{i}\left(k_{0}x-\omega t\right)},$$

(35)

immediately from (24) and (28), we obtain $k_1=k_0$ and

$$A_1=\frac{a_{0}g}{\mathrm{i}\omega }\frac{e^{-\mathrm{i}{k}_1(x_0+B)}}{cosh{k}_{1}h-\Omega k_{1}sinh\left(k_{1}h\right)},$$

(36)

where $a_0$ is the wave amplitude measured at some reference point $x=x_0$ and $k_0$ is the corresponding wavenumber. As for the remaining unknown coefficients associated with the solution for $\Phi$, namely $(B_n,C_m,D_m,E_n)$ in (28), (31), and (32), they are determined by the matching conditions typical for wave scattering problem by the method of eigenfunction expansions [4]. Specially, we shall require the continuity of both $\Phi$ and $\partial \Phi /\partial x$, i.e., pressure and velocity, across the vertical interfaces at $x=\pm B$. Unfortunately, the usual orthogonality of eigenfunctions that can be conveniently employed to deduce a linear system for unknown coefficients is no longer valid due to the presence of a muddy seabed, which can be readily understood by the nonzero vertical water particle velocity at the bottom of the water layer, as shown in (21). We reiterate that if $w_b=0$, the present model recovers the classical problem of wave scattering by a surface obstacle above a solid bed [4]. In the present study, to make the problem tractable, we shall adopt the idea of the so-called method of orthogonal mode coupling [11,24,25] as the treatment for the matching conditions. In other words, we define the following two inner products

$$\begin{array}{c}\langle f_n,f_m\rangle =\left\{\begin{array}{c}\left(1+{\Omega }^2{K}_n^2\right)\frac{sinh2k_{n}h}{4k_n}+\frac{\Omega \left(1-cosh2k_{n}h\right)}{2}+\frac{h\left(1-{\Omega }^2{k}_n^2\right)}{2},n=m\hfill \\ \begin{array}{cc}\hfill {\int }_{-h}^0{f}_n{f}_{m}dz-\frac{1+{\Omega }^2{k}_n{k}_m}{2}\frac{sinh\left(k_n+k_m\right)h}{k_n+k_m}-\Omega \left(1-cosh{k}_{m}hcosh{k}_{n}h\right)& \\ \hfill -\frac{1-{\Omega }^2{k}_n{k}_m}{2}\frac{sinh\left(k_n-k_m\right)h}{k_n-k_m},n\ne m\end{array}\hfill \end{array}\right.,\hfill \end{array}$$

(37)

$$\begin{array}{c}\langle p_n,p_m\rangle =\left\{\begin{array}{c}\frac{H}{2}+\frac{sin2K_{n}H}{4K_n},n=m\hfill \\ {\int }_{-h}^{-w}{p}_n{p}_{m}dz-\frac{K_{n}sin{K}_{n}Hcos{K}_{m}H-K_{m}cos{K}_{n}Hsin{K}_{m}H}{K_n^2-K_m^2},n\ne m\hfill \end{array}\right.,\hfill \end{array}$$

(38)

where $f_n=cosh{k}_n(z+h)-\Omega k_{n}sinh{k}_n(z+h)$ and $p_m=cos{K}_m(z+w)$ are those vertical functions that appear in the solution forms for $\Phi$, i.e., (28), (31), and (32). We note that $f_n$ and $p_m$ are orthogonal with respect to the above inner products [25],

$$\langle f_n,f_m\rangle =\langle p_n,p_m\rangle =0,n\ne m.$$

(39)

Using the above inner products, we obtain from the continuity of pressure at $x=-B$ and $x=B$

$$\begin{array}{cc}\hfill & \sum _{n=1}^N{B}_n{\Gamma }_{n,\tilde{m}}-\sum _{m=1}^M\left(C_m-D_m\right){\Delta }_{m,\tilde{m}}=-A_1{\Gamma }_{1,\tilde{m}},\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill & \sum _{n=1}^N{E}_n{\Gamma }_{n,\tilde{m}}-\sum _{m=1}^M\left(C_m+D_m\right){\Delta }_{m,\tilde{m}}=0,\hfill \end{array}$$

(41)

and from the matching of velocity

$$\begin{array}{cc}\hfill & \sum _{n=1}^N\mathrm{i}{k}_n{B}_n{\Theta }_{n,\tilde{n}}-\sum _{m=1}^M\left(C_m{K}_{m}tanh{K}_{m}B-D_m{K}_{m}coth{K}_{m}B\right){\Gamma }_{\tilde{n},m}=\mathrm{i}{k}_1{A}_1{\Gamma }_{1,\tilde{n}},\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & \sum _{n=1}^N\mathrm{i}{k}_n{E}_n{\Theta }_{n,\tilde{n}}-\sum _{m=1}^M\left(C_m{K}_{m}tanh{K}_{m}B+D_m{K}_{m}coth{K}_{m}B\right){\Gamma }_{\tilde{n},m,}=0,\hfill \end{array}$$

(43)

where

$$\begin{array}{cc}\hfill & {\Gamma }_{nm}=\frac{k_{n}sinh{k}_{n}H+K_{m}sin{K}_{m}H-\Omega k_n^2\left(cosh{k}_{n}H-cos{K}_{m}H\right)}{k_n^2+K_m^2},\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & {\Delta }_{nm}=\frac{K_{n}sin{K}_{n}Hcos{K}_{m}H-K_{m}cos{K}_{n}Hsin{K}_{m}H}{K_n^2-K_m^2},\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill & \begin{array}{cc}\hfill {\Theta }_{nm}=& \frac{1+{\Omega }^2{k}_n{k}_m}{2}\frac{sinh\left(k_n+k_m\right)h}{k_n+k_m}-\Omega \left(1-cosh{k}_{m}hcosh{k}_{n}h\right)\hfill \\ \hfill & -\frac{1-{\Omega }^2{k}_n{k}_m}{2}\frac{sinh\left(k_n-k_m\right)h}{k_n-k_m}.\hfill \end{array}\hfill \end{array}$$

(46)

In addition, $\tilde{m}=1,2,\cdots ,M$ and $\tilde{n}=1,2,\cdots ,N$. We note that the summations of infinite series are truncated at finite N and M terms with $N/M\approx h/H$ [4]. Hence, (40)–(43) now form a system of $(2M+2N)$ linear equations that can be solved numerically to obtain a total of $(2N+2M)$ unknown coefficients $(B_n,C_m,D_m,E_n)$.

4. Numerical Examples

We reiterate that by setting $d=0$, our model recovers the classical solution for wave scattering by a surface obstacle above a solid bed [4]. If there is no obstacle, i.e., $B=W=0$, the model reduces to the study of water waves over a layer of fluid mud [17]. We have checked our results with those reported in the literature [4,17].

With the surface obstacle and the muddy seabed both considered, to make a quick check, we compare our model predictions with the available numerical results [26] obtained using OpenFoam, which is a popular open-source CFD package. Figure 4 shows the wave amplitudes for incident, reflected, and transmitted waves. Visually, our results agree with the numerical solutions. The wave damping caused by the muddy seabed is evident. In Figure 5, we also compare the records of dynamic wave force acting on the bottom of the obstacle. The wave force, F, is obtained by integrating the dynamic pressure, which is calculated by the Bernoulli Equation (27), over the obstacle bottom face. The agreement between the model predictions and the numerical results is again reasonable, as can be seen in Figure 5. We notice that the dynamic force decreases slightly in time.

To understand the dynamic aspect of the impacts caused by the muddy seabed, in Figure 6, we show the distribution of dynamic pressure on the obstacle bottom face. Various values of wave frequency, ranging from long to short waves, are considered with other parameters unchanged. We observe from the figure that the pressure difference across the length of the obstacle is significant. The change in magnitude becomes more dramatic for longer waves. Although we consider only a fixed obstacle, the nonuniform distribution of dynamic pressure may cause a potential rotational failure of the structure, which requires more attention. In Figure 7, we consider the effects of obstacle submerged depth and the thickness of the mud layer. As can be seen, at a fixed value of $W/h$, the thicker the mud layer, the larger the difference in dynamic wave force, $\Delta F$. We note that the dimensionless wave force difference is defined as $\Delta F=(F_{\mathrm{mud}}-F_{\mathrm{solid}})/F_{\mathrm{solid}}$, showing the effects of a muddy seabed. On the other hand, when d is fixed, the change in wave force becomes larger as $W/h$ increases. We note that the passage below the obstacle becomes smaller for larger $W/h$, i.e., shallower depth below the obstacle. Hence, the cases with smaller H are more affected by the muddy seabed.

Finally, in Figure 8, we show the amplitudes of interfacial waves at the water–mud interface for short to long waves under several different values of obstacle submerged depth. When $W/h$ is fixed, the amplitude is larger for longer waves, which is in agreement with the common understanding that bottom conditions are more significant for long waves. On the other hand, for a given wave frequency, there exists a peak amplitude with respect to the relative obstacle submerged depth, $W/h$. Furthermore, the peak amplitude occurs at larger values of $W/h$ for shorter waves, which again agrees with the fact that long waves are more affected by the seabed.

5. Concluding Remarks

We present an analytical solution for wave scattering by a surface obstacle above a muddy seabed. The bottom mud is idealized as a Newtonian fluid, and the thickness of the mud layer is assumed to be smaller than the water depth. The present solution takes into account the combined effects of obstacle and fluid mud bottom. It reduces to the classical solutions reported in the literature if the mud layer or the obstacle is absent. Model predicted wave amplitude and wave force agree reasonably with the available numerical results. Effects of the mud layer thickness, the obstacle submerged depth, and the frequency of incident waves are discussed. We reiterate that the present model is limited to linear waves and the muddy seabed is considered thin. Although the fluid mud is assumed to act as a Newtonian fluid only, the model can also be applied to the case of a viscoelastic seabed. This is granted by the fact that for a linear problem, we only need to introduce a complex viscosity, ${\nu }_{me}={\nu }_m+\mathrm{i}\frac{G_m}{\omega {\rho }_m}$ with $G_m$ being the shear modulus of elasticity, and simply replace ${\nu }_m$ by ${\nu }_{me}$ in all our formulations [15,17]. To account for more complex rheological behaviors of mud, the present model can potentially be extended to incorporate the idea of multiple mud layers [28,29], providing that in each single layer, the mud acts as a viscoelastic fluid [22]. The layering structure can also be introduced in the water body if the possible effects of water density stratification are to be considered [30]. Finally, to relax the constant depth assumption in the present model, the idea of using a staircase seabed to approximate the varying bathymetry [31] is a potential candidate.