Numerical Investigation on Solitary Wave Interaction with a Vertical Cylinder over a Viscous Mud Bed

Abstract

This study investigated the hydrodynamics of a solitary wave passing a vertical cylinder over a viscous mud bed for the first time. A highly viscous Newtonian fluid was assumed as a simplified model for fluid mud. A three-dimensional numerical wave flume consisting of a fixed cylindrical structure and three viscous fluids—air, water, and mud—was constructed and validated. Numerical experiments were performed to investigate solitary wave interaction with a vertical cylinder over a viscous mud bed. Numerical results showed the mud surface deformation to be one order of magnitude smaller than the water surface deformation and their behaviors to be different: mud surface depressions occurred on the upstream and downstream sides of the cylinder, whereas mud surface elevations occurred on the lateral sides of the cylinder. This solitary wave induced scour pattern on a muddy seabed is different from that commonly observed on a sandy seabed. Water flow reversal near the water–mud interface was made more evident by the mud bed. Although the mud bed attenuated water waves, it nevertheless increased the total horizontal force and toppling moment exerted on the cylinder due to the wave-induced mud flow. These findings may be valuable to the design of marine structures on a muddy seabed and worthy of further investigation.

1. Introduction

The ocean makes up at least two-thirds of the Earth’s surface. As demand for energy grows, interest in harvesting resources in the ocean follow. Nowadays, abundant research on related topics can be found (see, e.g., a series of studies introduced by Correia et al. [1]). The design and construction of marine structures remain an important and relevant topic in ocean engineering. For example, scour around the foundations of marine structures, though a classical topic, still receives tremendous research efforts due to its complexity and practical value (see, e.g., a recent review by Fazeres-Ferradosa et al. [2]).

The placement of marine structures ranges from nearshore to offshore. In a nearshore environment, water waves are an important hydrodynamic factor to consider. Long waves (or shallow water waves), whose wavelengths are significantly larger than the water depth, are commonly assumed to study long-period waves such tsunamis, infragravity waves, and typhoon-induced swells in a nearshore environment.

Due to its clear mathematical definition and easiness to replicate, solitary waves are one of the most widely used long-wave models. For example, Grimshaw [3] derived a solitary wave solution in water of variable depth. Solitary waves also see numerous applications in studies on tsunamis (transient long waves), impacts of extreme waves on marine structures, nearshore sediment transport, swash zone dynamics, etc. (see e.g., Isaacson [4], Wang et al. [5], Yates and Wang [6], Tonkin et al. [7], Mo et al. [8], Lo and Liu [9], Pujara et al. [10], Larsen et al. [11,12], and Wu et al. [13]). The leading-order solution for the free surface displacement $\eta$ of a solitary wave reads

$$\eta \left(t\right)=H{\mathrm{sech}}^2\left[\sqrt{\frac{3H}{4h^3}}(x-ct)\right],$$

(1)

in which H is the wave height, h is the still water depth, g is the gravitational acceleration, t is time, x is the horizontal spatial coordinate, and

$$c=\sqrt{g(H+h)}$$

(2)

is the wave speed. Strictly speaking, solitary waves have infinitely long wavelength. However, in practice, an effective wavelength l is often used:

$$l=4\pi \sqrt{\frac{h^3}{3H}}.$$

(3)

To minimize wave drag in all directions, the foundation piles of marine structures usually have circular cross sections (see, e.g., Jamalabadi and Oveisi [14]). In studies on wave and current interactions with marine structures, vertical cylinders are often assumed as an approximate shape of foundation piles. For example, Isaacson [4] derived the analytical solution for solitary wave diffraction around a large cylinder; Wang et al. [5] performed numerical simulations of the same problem based on a depth-averaged wave model; Yates and Wang [6] conducted laboratory experiments on the same problem and compared the wave runup around the cylinder and the forces exerted on the cylinder with numerical results; Mo et al. [8] performed three-dimensional (3D) numerical simulations of wave–cylinder interaction for both regular waves and solitary waves, and validated the numerical results through large-scale laboratory experiments; Tonkin et al. [7] experimentally investigated tsunami-induced scour around a coastal cylindrical structure; Larsen et al. [11] numerically investigated tsunami-induced scour around a monopile structure, representative of offshore wind turbine foundations in moderate water depths; Larsen et al. [12] performed laboratory experiments on the same problem and validated the numerical results in Larsen et al. [11]; Li et al. [15] conducted 3D numerical simulations to mimic the Tonkin et al. [7] laboratory experiments; Li et al. [16] conducted 3D numerical simulations to study tsunami-induced scour around a bridge pier.

In modeling, seabeds can be loosely divided into three categories: rigid, sandy, and muddy. Conventionally, either a rigid seabed or a sandy seabed is assumed in studies on wave-structure interaction, such as the above-mentioned examples. Simply put, a rigid seabed consists of bedrock and boulders, often modeled as an impermeable solid bottom boundary. A sandy seabed consists of granular (large-grain) materials like sand and gravel, often modeled as a porous media whose pore pressure effects must be considered. As pointed out by Healy et al. [17], muddy seabeds exist in many parts of the world. While the scour of sandy seabeds has been studied extensively (see, e.g., Sumer et al. [18]), studies on the scour of muddy seabeds are scarce. A muddy seabed consists of cohesive (fine-grain) sediments such as mud and clay, whose inter-particle adhesive forces must be considered and a rheological model is often assumed. Studies have indicated that cohesive sediments demonstrate different scour and erosion characteristics from noncohesive granular materials and warrant further investigation (see, e.g., Schindler et al. [19], Harris and Whitehouse [20]). While scour protection countermeasures in uniform noncohesive soils have been widely studied (see e.g., Fazeres-Ferradosa et al. [21]), knowledge on the behavior of scour and scour protections in cohesive soils or layered noncohesive soils is lacking (see, e.g., Fazeres-Ferradosa et al. [2], Porter et al. [22]). This study seeks to improve the understanding of wave-structure interaction over a muddy seabed (i.e., cohesive sediments).

Wave–mud interaction, especially the attenuation of waves by a muddy seabed, has been studied by various researchers. As the most simplified model, the mud is commonly assumed to be a highly viscous Newtonian fluid mud. Gade [23] derived analytical solutions based on a two-fluid model, in which the upper layer is a shallow-depth inviscid fluid (water) and the lower layer is a Newtonian viscous fluid mud, and compared the results with experimental measurements. Dalrymple and Liu [24] extended the Gade [23] model to any water depth and considered the viscosity of water, and found that high wave attenuation rates are possible when the mud layer is thick or when the mud layer is about as thick as its boundary layer. By considering the water boundary layer near the water–mud interface, the mud boundary layer near the water–mud interface, and mud boundary layer near the rigid bottom, Jiang and Zhao [25] analyzed the attenuation of solitary waves propagating over a Newtonian fluid-mud seabed, and conducted laboratory experiments to validate the analytical results. Based on the Boussinesq wave model, Liu and Chan [26] derived a set of depth-integrated equations for long-wave propagation over a Newtonian fluid-mud seabed, and calculated the damping rates for periodic and solitary waves. Park et al. [27] conducted laboratory experiments to validate the analytical solutions derived by Liu and Chan [26].

Although the high-viscosity Newtonian fluid mud model is often assumed in existing studies, it may not be the most accurate for real cohesive sediments like clay and mud. Measurements have shown real cohesive sediments to demonstrate shear thinning characteristics, where the shear stress varies with the shear rate (see e.g., Hsu et al. [28], Guo et al. [29,30]). Therefore, non-Newtonian fluid mud models may be more fitting, although also more complex, options. As one of the simplest non-Newtonian fluid models, and a direct extension of the Newtonian fluid model, the Bingham-plastic model is a commonly assumed non-Newtonian model for mud; see, e.g., Mei and Liu [31], Chan and Liu [32], Hsu et al. [28], and Guo et al. [33]. Other rheological models are also possible: Mei et al. [34] assumed a viscoelastic model for mud, and Oveisy et al. [35] assumed a visco-elastic-plastic model for mud.

To better understand the factors contributing to wave attenuation, the viscous damping of a solitary wave due to its boundary layer near a rigid bed has also attracted attention. Keulegan [36] first derived a leading-order model, Liu and Orfila [37] extended the Keulegan [36] model to higher order, and Liu et al. [38] extended the Liu and Orfila [37] model to include nonlinear terms and conducted laboratory experiments to validate the analytical findings. Flow reversal, i.e., local velocity pointing in the opposite direction to the freestream velocity, was observed in the boundary layer of a solitary wave over a rigid bed.

The literature review above reveals abundant studies on wave–mud interaction as well as wave–cylinder interaction on either a rigid seabed or a sandy seabed. However, to the best of the authors’ knowledge, no studies on wave–mud–cylinder interaction have been published. In designing monopile foundations in areas known to have a muddy seabed, simultaneous consideration of waves, cylindrical structure, and mud may better resemble site conditions. For example, the proposed offshore wind farms in the Taiwan Strait are mostly located in water depths of 50 m or less. Long-period typhoon-induced swells are not uncommon in this region. In addition, studies have found the marine sediments in the Taiwan Strait to consist of varied mixtures of coarse-grain sediments like sand and fine-grain sediments like clay (see e.g., Xu et al. [39], Huh et al. [40]). Therefore, investigations on long wave interaction with a vertical cylinder over a fluid mud bed shall be of both practical value and scientific value.

This study seeks to perform a numerical investigation on wave–mud–cylinder interaction. As a preliminary investigation of a previously unexamined problem, this study adopts the widely used solitary wave as the incident long wave. A highly viscous Newtonian fluid is assumed as a simplified model for fluid mud, using the material properties of the “fluid mud” (a silicone fluid) from the Park et al. [27] experiments. The 3D multiphase flow problem consisting of three incompressible Newtonian fluids—air, water, and fluid mud—is solved using the opens-source CFD (computational fluid dynamics) package OpenFOAM. The basic Newtonian fluid mud model is assumed in this study because it is the only model for which both analytical studies and detailed experimental data on wave–mud interaction are readily available for checking the simulated results. In order to avoid extrapolating the numerical simulations beyond the validated range of applicability, presently the numerical experiments are scaled based on relevant benchmark laboratory experiments (i.e., Park et al. [27] and Yates and Wang [6]) rather than real site conditions. Extending the numerical configuration to consider non-Newtonian fluid mud models and real site conditions is beyond the scope of this study.

This paper is organized as follows: in Section 2, the OpenFOAM configuration used in this study is outlined; in Section 3, the numerical model is validated using two relevant laboratory experiments; in Section 4, the design of the numerical experiments and mesh convergence tests are presented; in Section 5, the numerical results are discussed, including the water surface displacement, the mud surface displacement, the velocity field, and the forces and moments exerted on the cylinder; lastly, in Section 6, we conclude this study and remark on potential future research topics.

2. Numerical Model

To construct a three-dimensional numerical wave flume including water, viscous fluid mud, a vertical cylinder, and an incident solitary wave, this study utilizes the multiphase flow solver “multiphaseInterFoam” of the open-source CFD package OpenFOAM (Weller et al. [41]). The wavemaking boundary condition from the olaFlow module (Higuera [42]) is used to input a solitary wave into the computational domain. OpenFOAM is based on the finite volume method (Greenshields [43]), which calculates the cell-averaged values of each physical variable. The Euler implicit scheme is adopted to advance in time; the OpenFOAM built-in “PIMPLE algorithm” is adopted to solve the velocity-pressure coupling problem for incompressible flows. Since the present investigation does not consider breaking waves, no turbulence model is used.

Air–water and water–mud interfaces are modeled using the Volume of Fluid (VOF, Hirt and Nichols [44]) method. The VOF method assumes the fluids to be immiscible with each other, and assigns a coefficient ranging from 0 to 1 to denote the volume fraction of a fluid in each numerical cell. In this study, the water surface is defined as the linearly interpolated vertical elevation where the volume fraction of water is $0.5$; the mud surface is defined as the linearly interpolated vertical elevation where the volume fraction of mud is $0.5$.

OpenFOAM solves for the cell-averaged flow velocities and pressure, and outputs the results at each cell center. To estimate the wave forces acting on a vertical cylinder, we integrate pressure over the surface area of the cylinder. For example, to estimate the horizontal wave force acting on a vertical cylinder in the x-direction (i.e., the streamwise direction), we assume that pressure is constant within each cell. For each cell in contact with the surface of the cylinder, we multiply the cell pressure with its effective surface area projected in the x-direction. Then, we sum up the force contribution of each cell to obtain the net horizontal wave force acting on the cylinder.

The main computational domain of this study is a three-dimensional wave flume. As sketched in Figure 1, a slot is specified on the flume bottom in the middle of the flume, to be filled with viscous fluid mud. In addition, a vertical cylinder is fixed at the center of the flume. Of course, depending on the problem to be simulated, the exact configuration of the numerical wave flume varies. The no slip condition is imposed on all solid boundaries: “wall”, “bottom”, “front”, “back”, and “cylinder”, as labeled in Figure 1; an incident solitary wave is specified on the “inlet” boundary; the “atmosphere” boundary is open to the atmosphere. Since the right wall of the flume is a solid boundary, wave reflection occurs after the incident wave reaches the right wall. Nonetheless, we analyze and present only the results before the arrival of the reflected waves. Therefore, wave reflection from the right wall is irrelevant in this study.

3. Validation Tests

OpenFOAM in combination with the olaFlow module has been widely applied to study free surface flows (see, e.g., the references provided by Higuera [42]), including the generation and propagation of solitary waves. However, to the best of our knowledge, it has not been applied to investigate wave–structure interaction over a muddy seabed. Hence, we find it necessary to perform and present validation tests relevant to the problem of interest—solitary wave interaction with a vertical cylinder over a viscous mud bed. Unfortunately, no experimental measurements resembling this configuration are available for a direct validation test. As second-best options, we select two validation tests, each resembling part of the wave–cylinder–mud problem. The two tests are: a solitary wave propagating over a viscous mud bed, and a solitary wave passing a vertical cylinder over a rigid bottom.

3.1. Solitary Wave Propagating over a Viscous Mud Bed

This section presents the validation test for a solitary wave propagating over a viscous mud bed, based on the Park et al. [27] laboratory experiment, which was conducted in a physical wave flume of length of 32 m, width $0.6$ m, and depth $0.9$ m. A piston-type wavemaker was installed on the left side of the flume. A slot of length $0.97$ m, width $0.6$ m, and height $0.017$ m was placed at a distance of $12.3$ m away from the left side of the flume. To mimic a muddy seabed, the slot was filled with a highly viscous Newtonian silicone fluid. The silicone fluid has a kinematic viscosity of 5.24 × 10⁻³ m²s⁻¹ and a density of 1050 kg/m³. It is transparent, allowing for optical measurements of the flow. Outside of the mud slot, the flume was filled with water to a depth of $0.1$ m. A solitary wave with a wave height of $0.019$ m and an effective wavelength of $1.72$ m was generated by the wavemaker as the incident wave. According to Park et al. [27], based on the Reynolds number similarity, this configuration approximately translates to a prototype water depth of 10 m over a viscous mud with a kinematic viscosity of 5.36 m²s⁻¹. In comparison, Gade [23] estimated the viscosity of the top-layer sediments in the Mississippi River Delta, USA to be between 0.1 m²s⁻¹ and 1 m²s⁻¹.

In our numerical simulations, the kinematic viscosity of water was set to be 10⁻⁶ m²s⁻¹, and the density of water was set to be 1000 kg/m³. To increase computational efficiency, a shorter wave flume was used. As sketched in Figure 2, the numerical wave flume is 8.97 m long. The origin of the coordinate system is defined at the bottom left of the numerical wave flume; i.e., x = 0 m is defined at the left side of the flume, and z = 0 m is defined at the bottom of the left side of the flume.

Park et al. [27] measured the water surface displacement, the water flow velocities, the mud surface displacement, and the mud flow velocities. In Figure 3, we compare the numerical results at $x=4.43$ m for the water surface and the horizontal water flow velocity (at $z=0.012$ m) with experimental measurements. Numerical results based on three different mesh sizes are shown: coarse mesh ($\Delta x=1$ cm, $\Delta y=0.5$ cm), medium mesh ($\Delta x=0.5$ cm, $\Delta y=0.25$ cm), and fine mesh ($\Delta x=0.25$ cm, $\Delta y=0.125$ cm). The Grimshaw [3] solitary wave analytical solutions are also plotted. The same symbols and nondimensionalization factors as in Park et al. [27] are used. Thus, for brevity they need not be reiterated here. It can be seen that the results agree well with each other. The coarse-mesh results slightly underestimate the peak values, whereas the medium-mesh results are nearly identical to the fine-mesh results, indicating that the medium mesh is sufficient for simulating the water wave in this problem.

In Figure 4, we compare the numerical results at $x=4.43$ m for the horizontal mud flow velocity (at $z=-0.01$ m), and the mud flow velocity profiles with experimental measurements. The analytical solutions from Liu and Chan [26] are also plotted. Again, the same symbols and nondimensionalization factors as in Park et al. [27] are used, which need not be reiterated here. In Figure 4a, numerical results based on three different mesh sizes are shown for the horizontal mud flow velocity. Although the medium mesh is sufficient for simulating the water wave, to accurately simulate the mud flow is more demanding—Figure 4a shows that the fine mesh is needed for the numerical results to agree with experimental measurements. Therefore, in Figure 4b, only the fine-mesh numerical results are shown.

Overall, good agreement between the fine-mesh numerical results and the experimental measurements is achieved. Therefore, we have validated the capability of the numerical model to simulate a solitary wave propagating over a viscous mud bed.

3.2. Solitary Wave Passing a Vertical Cylinder over a Rigid Bottom

This section presents the validation test for a solitary wave passing a fixed vertical cylinder over a rigid bottom, based on the Yates and Wang [6] laboratory experiment, which was conducted in a physical wave flume of length $7.6$ m, width $0.76$ m, and height $0.6$ m. A piston-type wavemaker was installed on the left side of the flume. Along the center line and in the middle of the flume, a vertical cylinder of radius $0.0635$ m and height $0.6$ m was installed. The still water depth was $0.04$ m, and the incident wave height of the solitary wave was $0.016$ m. According to Yates and Wang [6], this configuration falls in the large-body and shallow-water regime. The diffraction parameter is

$$\frac{D}{l}=0.277,$$

(4)

in which D is the cylinder diameter, and l is the effective wavelength of a solitary wave, Equation (3). Based on the criterion derived by Isaacson [4], solitary wave diffraction by the cylinder is important in this experiment. For a solitary wave, Yates and Wang [6] calculated the Keulegan–Carpenter number to be

$$K_C=\frac{H}{h}\frac{l}{D}=1.445,$$

(5)

in which H is the incident wave height and h is the still water depth as defined in (1). For $K_C<2$, drag forces can be neglected in comparison to inertia forces (see, e.g., Sarpkaya and Isaacson [45] and Yates and Wang [6]).

The three-dimensional numerical wave flume used for this validation test is sketched in Figure 5. The origin of the coordinate system is defined at the center of the base of the cylinder. x points in the streamwise direction, y points in the lateral direction, and z points in the vertical direction. To increase computational efficiency, a shorter numerical wave flume of length $3.8$ m was used.

Yates and Wang [6] measured the water surface displacement and the horizontal wave force exerted on the cylinder. The numerical results are compared with experimental measurements for the water surface displacement at $(x,y)=(-0.066,0)$ m in Figure 6a, and for the horizontal wave force exerted on the cylinder in Figure 6b. The same symbols and nondimensionalization factors as in Yates and Wang [6] are used. Hence, for brevity they need not be reiterated here. Numerical results based on three different mesh sizes are shown: coarse mesh (smallest cell size: $\Delta x=0.5$ cm, $\Delta y=0.8$ cm, $\Delta z=0.2$ cm), medium mesh (smallest cell size: $\Delta x=0.375$ cm, $\Delta y=0.6$ cm, $\Delta z=0.15$ cm), and fine mesh (smallest cell size: $\Delta x=0.25$ cm, $\Delta y=0.4$ cm, $\Delta z=0.1$ cm). For these three mesh sizes, the numerical results are not very sensitive to the mesh sizing, as can be seen in Figure 6.

The water surface results, i.e., Figure 6a, shows a slight difference in the wave shape between the numerical results and the experimental measurements. This is likely due to the less-than-ideal solitary wave generation in the laboratory, since this deviation from experimental measurements also shows in the numerical simulation of this problem performed by Wang et al. [5]. Aside from the slight shape difference, the overall agreement between the numerical results and the experimental measurements is good. Therefore, we have validated the capability of the numerical model to simulate a solitary wave passing a fixed vertical cylinder over a rigid bottom.

4. Setup of Numerical Investigation

Having validated the numerical model for simulating wave–mud interaction and wave–cylinder interaction, we shall begin the numerical investigation on the wave–mud–cylinder interaction. This section is dedicated to describing the numerical configuration. The results will be presented in the next section.

4.1. Numerical Experiments

Adopting the concept of control variables, we devised four numerical experiments. The configurations of these four numerical experiments are summarized in

Table 1. Configurations of the four numerical tests.

Name	Test1	Test2	Test3	Test4
Incident wave	solitary wave ($H/h=0.4$)	solitary wave ($H/h=0.4$)	solitary wave ($H/h=0.4$)	solitary wave ($H/h=0.4$)
Cylinder	yes	yes	yes	no
Slot	yes	yes	no	yes
Material in slot	highly viscous Newtonian fluid	water	N/A	highly viscous Newtonian fluid

. Comparing the results from these four numerical experiments shall reveal the influences of the cylinder, of the bottom slot, and of the fill material placed in the bottom slot. The same incident wave, a solitary wave with $H/h=0.4$, was used in all four numerical experiments.

As an example, the numerical wave flume used for test1 is sketched in Figure 7. The origin of the coordinate system is defined at the center of the cylinder and on the initial water–mud interface. Therefore, in still water, water occupies $z\ge 0$, and the viscous fluid mud occupies $z<0$. To ensure that the numerical simulations are within the validated range of applicability, scaling and parameter choice in the present investigation are based purely on existing laboratory experiments.

As in the validation test based on the Yates and Wang [6] experiments, the still water depth is $0.04$ m, the incident solitary wave has a wave height of $0.016$ m, and the effective wavelength is $0.4589$ m. Since these scales are the same as those in Yates and Wang [6], the wave–cylinder interaction falls in the large-body and shallow-water regime. The diffraction parameter is $0.277$, Equation (4)—wave diffraction by the cylinder is significant; the Keulegan–Carpenter number is $1.445$, Equation (5)—drag forces are negligible in comparison to inertia forces.

As indicated in Figure 7, a bottom slot of length $0.4$ m, width $0.76$ m, and depth $0.02$ m is installed in the middle of the flume. The vertical cylinder extends all the way to the bottom of the slot; i.e., $z=-0.02$ m. The slot is filled with the highly viscous silicone fluid used by Park et al. [27]. Its kinematic viscosity is 5.24 × 10⁻³ m²s⁻¹, and its density is 1050 kg/m³. The kinematic viscosity of water is set to be 10⁻⁶ m²s⁻¹, and the density of water is set to be 1000 kg/m³. Since the water depth and the mud slot depth are both different from the Park et al. [27] experiment, the prototype scales are also different. Based on the Reynolds number similarity proposed by Park et al. [27], our numerical configuration approximately translates to a prototype water depth of 10 m over a viscous mud with a kinematic viscosity of 20.71 m²s⁻¹, about four times more viscous than the prototype viscosity of the Park et al. [27] experiment.

4.2. Mesh Configuration

As sketched in Figure 8, a structured yet nonuniform mesh is used. Since our main interests are the fluid interfaces and flows near the slot and the cylinder, a graded mesh is utilized. Mesh spacing in the x-direction is coarser near the two ends of the flume and gets finer near the slot; mesh spacing in the z-direction is finer near the air–water interface and the water–mud interface; mesh spacing in the y-direction is uniform. As can be seen in Figure 8a, the mesh is curvilinear near the cylinder; elsewhere, the mesh is rectilinear.

Mesh convergence tests for test1 were conducted based on this mesh configuration. In the coarse mesh, the smallest cell size is $(\Delta x,\Delta y,\Delta z)=(0.5,0.8,0.2)$ cm. In the medium mesh, the smallest cell size is $(\Delta x,\Delta y,\Delta z)=(0.375,0.6,0.15)$ cm. In the fine mesh, the smallest cell size is $(\Delta x,\Delta y,\Delta z)=(0.25,0.4,0.1)$ cm. Comparisons of the numerical results for the water surface displacement at $(x,y)=(-0.1325,0)$ m, the mud surface displacement at $(x,y)=(-0.1325,0)$ m, the velocity profile at $(x,y)=(-0.1325,0)$ m, and the horizontal force exerted on the cylinder are shown in Figure 9. While the results based on differently sized meshes agree well with each other for the water surface displacement and the horizontal force, discrepancies are noticeable for the mud surface displacement and the velocity profile. We note that the mud surface displacement is one order of magnitude smaller than the water surface displacement. In addition, the mud surface displacement is smaller than the vertical mesh spacing, and the VOF method is not designed to precisely track the interfaces. Despite the discrepancies, the trends and characteristics of the results agree with each other regardless of the mesh size.

Numerical experiments with an event duration of 5 s were simulated on the Intel i7-9800X CPU ($3.80$ GHz $\times 8$). The coarse-mesh computation time is approximately 10 h, the medium-mesh computation time is approximately 7 days, and the fine-mesh computation time is approximately 54 days. As this study is a preliminary investigation focusing more on the qualitative understanding, rather than precise quantitative measurements, of an underexplored problem, the coarse mesh was adopted due to its feasibility and computational efficiency. More details on the mesh configuration based on the coarse mesh are tabulated in

Table 2. Descriptions of the numerical meshes used in the four numerical experiments.

Name	Test1	Test2	Test3	Test4
Flume size	3.80 m long, 0.76 m wide, 0.08 m high	same as test1	same as test1	same as test1
Slot size	0.40 m long, 0.76 m wide, 0.02 m high	same as test1	N/A	same as test1
Cylinder	radius 0.0635 m	same as test1	same as test1	N/A
Min. mesh	$(\Delta x,\Delta y,\Delta z)=(0.5,0.8,0.2)$ cm	same as test1	same as test1	same as test1
$\Delta x$	$\Delta x=0.5$ cm for $-0.2\le x\le 0.2 \mathrm{m}$; increases to $\Delta x=0.8$ cm at $x=\pm$1.9 m	same as test1	same as test1 (w/o slot)	same as test1 (w/o cylinder)
$\Delta y$	$\Delta y=0.8$ cm for rectangular cells; $\Delta y\simeq 0.8$ cm for curvilinear cells	same as test1	same as test1 (w/o slot)	same as test1 (w/o cylinder)
$\Delta z$	$\Delta z=0.2$ cm near the air–water and water–mud interfaces; increases to $\Delta z=0.5$ cm away from the interfaces	same as test1	same as test1 (w/o slot)	same as test1 (w/o cylinder)
Mesh near cylinder	4 concentric layers from the cylinder; each layer consists of 76 curvilinear cells of roughly equal size as nearby rectangular cells	same as test1	same as test1 (w/o slot)	same as test1 (w/o cylinder)
No. of cells	1,609,884	same as test1	1,557,468	1,628,110

for the four numerical experiments to be investigated.

5. Results

Having described the mesh configuration, we now present and discuss the numerical results based on the coarse mesh for the four numerical experiments listed in Table 1. The incident solitary wave travels from negative x towards positive x. $t=0$ s is defined as the instant when the freestream (i.e., away from the cylinder) wave crest passes $x=0$. Therefore, the wave is upstream of the cylinder for $t<0$ s, and downstream of the cylinder for $t>0$ s.

5.1. Water Surface and Mud Surface

Snapshots of the water surface displacement and the mud surface displacement are shown in Figure 10 from $t=-0.2$ s to $t=0.4$ s. In Figure 10a at $t=-0.2$ s, the wave crest approaches the cylinder. Both the water surface and the mud surface rise upstream of the cylinder. In Figure 10b at $t=0$ s, the wave crest reaches the lateral axis of the cylinder ($x=0$ m). Water surface upstream of the cylinder and on the lateral sides of the cylinder rises. The mud surface rises on the lateral sides of the cylinder. However, depressions occur on the mud surface upstream of the cylinder. In Figure 10c at $t=0.2$ s, the wave is just downstream of the cylinder. Wave runup on the cylinder can be observed on the downstream side. Due to wave scattering, an area of water surface depressions shows on the upstream and lateral sides of the cylinder. The mud surface also shows a depression upstream of the cylinder. However, differently from the water surface, on the lateral sides, the mud surface rises; a slight mud surface depression also shows immediately downstream of the cylinder. In Figure 10d at $t=0.4$ s, the wave has just left the cylinder. The water surface depression due to wave scattering continues radiating outwards from the cylinder. The mud surface depression upstream of the cylinder remains, the mud surface elevation on the lateral sides of the cylinder rises near to a maximum, and the mud surface depression downstream of the cylinder also grows near to a maximum extent.

From all panels of Figure 10, it can be seen that the mud surface displacement is one order of magnitude smaller than the water surface displacement; i.e., $O\left(1\right)$ mm vs. $O\left(1\right)$ cm. The water wave height is slightly reduced by the mud bed as expected per the derivation by Liu and Chan [26]. Nonetheless, the overall solitary wave scattering pattern due to a vertical cylinder remains similar to that predicted by Wang et al. [5]. Although the mud surface displacement is induced by the water waves, it deforms differently from the water surface displacement—on the upstream side of the cylinder, the mud surface first rises, and then depresses; on the lateral sides of the cylinder, the mud surface rises; on the downstream side of the cylinder, the surface depresses. These seabed deformation characteristics differ significantly from typical solitary-like wave scour on a sandy seabed around a cylinder, which tend to show a scour hole around the cylinder, and sediment deposits some distance downstream of the cylinder (see e.g., Larsen et al. [12], Li et al. [16]). Albeit the present study is a preliminary investigation assuming a highly idealized fluid mud model (i.e., a highly viscous Newtonian fluid), it nonetheless reveals the possibility of achieving a different scour pattern in the event of a muddy seabed, as opposed to the commonly assumed sandy seabed. As the scour problem is important in both foundation design and sediment transport, further investigations assuming a more realistic mud model may be needed.

5.2. Velocity Field

Just like the water surface, the freestream (i.e., away from the mud bed) water flow velocity is only slightly affected by the mud bed. Since the freestream water flow velocity is similar to that in existing studies on a solitary wave passing a vertical cylinder over a rigid bottom (e.g., Isaacson [4], Wang et al. [5]), our focus is on the water flow velocity near the mud bed and the mud flow velocity, which are understudied phenomena that warrant further discussions.

Snapshots of the horizontal water flow field at $z=0.003$ m (near the water–mud interface or the water–solid bottom boundary) and $t=0$ s (when the wave crest reaches $x=0$ m) are shown in Figure 11 for the four numerical experiments listed in Table 1. In Figure 11b,c, where no mud is present, a typical flow-around-a-cylinder pattern shows, with the horizontal velocities mostly pointing towards the positive direction. However, in Figure 11a,d, where a mud bed is present, flow reversal can be clearly observed: the horizontal velocities are zero near $x\simeq -0.06$ m, positive for $x\gtrsim -0.06$ m, and negative for $x\lesssim -0.06$ m. In addition, the flow speeds in Figure 11a,d are noticeably smaller than those in Figure 11b,c—namely, the water–mud interface slows down the flow more than the water–solid interface. On the other hand, due to an increase in the effective local water depth in test3, to satisfy mass conservation the flow speeds in Figure 11b appear slower than those in Figure 11c. By comparing Figure 11a with Figure 11d, it can also be seen that the existence of the cylinder increases the reverse flow speeds immediately upstream of the cylinder.

To show the vertical distributions of the horizontal velocities, the velocity profiles at three representative locations are examined—the upstream location $A(x=-0.1325,y=0)$ m, the freestream location $B(x=0,y=0.2240)$ m, and the downstream location $C(x=0.1325,y=0)$ m, as sketched in the bottom right of Figure 12. Velocity profiles from the four numerical experiments listed in Table 1 are plotted in Figure 12 at $t=0$ s.

At the upstream location A, both the test1 results and the test4 results show a reversal of water flow—the horizontal velocity changes sign near $z\simeq 0.006$ m in test1, and the horizontal velocity changes sign near $z\simeq 0.004$ m in test4. This difference in the elevation of water flow reversal is most likely due to the damming effects of the cylinder—the only configuration difference between test1 and test4. In addition, the mud flow velocity also points in the negative direction, and its overall magnitude is smaller than the water flow.

In the results for test2 and test3, where no mud is present, water flow reversal is not evident under the present mesh resolution, indicating that the presence of a mud bed enlarges the region in which flow reversal may occur.

At the freestream location $B(x=0,y=0.2240)$ m and the downstream location $C(x=0.1325,y=0)$ m, flow velocities for test1 are all positive. It can be seen that water flow is slowed down near the water–mud interface, and that the mud flow is noticeably slower than the water flow.

In test1, due to the presence of both a vertical cylinder and a viscous mud bed, a clockwise vortex on the $xz$-plane forms immediately upstream of the cylinder when the local water surface starts to fall. This phenomenon resembles the horseshoe vortex commonly observed in the scour problem. A snapshot at $t=0$ s is shown in Figure 13a for test1. As shown in Figure 13b for test3, this phenomenon is not evident in the absence of a mud bed.

After inspecting the velocity results in Figure 11, Figure 12 and Figure 13, it is found that the mud bed not only slows down the water flow near the water–mud interface, but also appears to facilitate water flow reversal. Hence, in the next section we shall examine the flow reversal phenomenon further.

5.3. Flow Reversal

Water flow reversal near the water–mud interface is observed in the numerical results. Liu et al. [38] pointed out that flow reversal tends to occur within boundary layers. Therefore, the study on water flow reversal is closely linked to the study of wave-induced boundary layer flow, on which various studies exist; e.g., Keulegan [36], Liu and Orfila [37], Liu et al. [38] examined the water boundary layer of a solitary wave over a rigid bottom, and Dalrymple and Liu [24], Jiang and Zhao [25], Liu and Chan [26], Park et al. [27] examined the mud boundary layer induced by a solitary wave traveling over a mud bed.

Existing studies did not consider the boundary layer induced by a solitary wave passing a vertical cylinder, nor was particular attention paid to the water flow reversal phenomenon for waves traveling over a muddy seabed. Hence, this section shall further discuss this phenomenon based on the newly obtained numerical results. We stress that this study is intended as a preliminary investigation of an underexplored problem, thus an elaborate boundary layer flow analysis is beyond the scope of this study.

To better visualize the location where flow reversal occurs, we define the vertical location where flow reversal occurs as the vertical location where the horizontal velocity in the x-direction changes sign; a linear velocity distribution between two vertically adjacent numerical cells is assumed. An example comparing the test1 results and the test4 results at $t=0.4$ s is shown in Figure 14. Because of the existence of a mud bed, water flow reversal is facilitated and captured in both cases using the present mesh resolution. The effects due to the cylinder are evident—it allows water flow reversal to occur further upstream near the cylinder and increases the vertical location of occurrence on the lateral sides of the cylinder. However, water flow reversal does not occur immediately downstream of the cylinder.

Next, we examine the simpler case test4 more closely. As suggested by Liu and Chan [26], in wave–mud interaction, a pressure gradient can be used to explain the behaviors of the mud flow and the water flow near the water–mud interface. Therefore, we shall utilize the numerical results for test4 to examine the relations between the water surface, the pressure gradient, and water flow reversal. A two-point central difference is adopted to estimate the horizontal pressure gradient $\partial p/\partial x$. Snapshots at $t=-0.35$ s and $t=0.05$ s are provided in Figure 15.

As can be seen in Figure 15, the pressure gradient is negative ahead of the wave crest and is positive behind the wave crest. Flow reversal occurs in the presence of a positive pressure gradient, which is referred to as the wave deceleration phase by Liu et al. [38]. It can also be seen that regardless of whether flow reversal occurs or not, water flow velocities near the water–mud interface are noticeably smaller than those in the freestream. The near-interface region (approximately $0\le z\le 3$ mm), in which water flow slows down, is significantly larger than the solitary wave boundary layer thickness of $0.87$ mm predicted by Liu and Chan [26]. Hence, the numerical results seem to imply that the mud bed increases the thickness of the water boundary layer induced by a solitary wave propagating over a viscous mud bed; in Liu and Chan [26], the boundary layer thickness is assumed to be the same regardless of whether the bottom is rigid or a viscous mud bed. Another possible explanation is that the water flow near the water–mud interface slows down due to the momentum transfer to the mud flow.

To summarize, the flow reversal phenomenon is further discussed in this section. Although the flow reversal under a solitary wave over a rigid bottom has been well studied by Liu and Orfila [37] and Liu et al. [38], that over a viscous mud bed has not been studied carefully. Liu and Chan [26] assumed the solitary wave boundary layer thickness on a viscous mud bed to be the same as that on a rigid bottom. However, newly obtained numerical results suggest that may be untrue—a vast difference between the $0.87$ mm thickness predicted by existing theory and the ∼3 mm thickness observed in the numerical results. As a preliminary investigation, the present mesh resolution in the vertical direction is at best 2 mm, which is inadequate for precisely describing boundary layer flow characteristics. Therefore, future investigations focusing more on this subject are needed.

5.4. Mud Flow

Per the expressions provided by Park et al. [27], the ratio of the mud layer thickness to the mud boundary layer thickness in the present numerical experiments is estimated to be $0.322$. The ratio was $0.178$ in the Park et al. [27] laboratory experiments. As the ratio is smaller than one, the mud flow is purely a boundary layer flow and has an approximately parabolic profile according to Park et al. [27]. As can be seen in Figure 12, this is indeed true.

Since the mud flow in the numerical experiments is a boundary layer flow, it is primarily driven by the viscous effects and the pressure gradient. As can be seen in Figure 15, the mud flow corresponds nearly instantly to the pressure gradient—the mud moves forward in the presence of a negative pressure gradient ahead of the wave crest, decelerates quickly, and eventually moves in the opposite direction in the presence of a positive pressure gradient behind the wave crest. In addition, both the mud surface and the mud flow velocities demonstrate oscillatory characteristics. While the water surface and the freestream velocities of a solitary wave are strictly nonnegative, those of the mud oscillate between positive and negative values.

5.5. Forces and Moments

Wave forces and moments exerted on a vertical cylinder are of practical engineering value. The horizontal forces exerted on the cylinder from the three numerical experiments listed in Table 1 that have a vertical cylinder are compared in Figure 16a; the toppling moments are compared in Figure 16b. A positive force means a net force in the positive x-direction; a negative force means a net force in the negative x-direction. Mud flow induced forces are included in the force calculation whenever applicable (i.e., in test1 and test2). Toppling moments are taken at the base of the cylinder along the pitch axis; i.e., along the $(x=0,y,z=-0.02)$ line for test1 and test2, and along the y-axis for test3, per the coordinate system sketched in Figure 7. A positive toppling moment means a moment that tends to topple the cylinder towards the downstream direction; a negative toppling moment means a moment that tends to topple the cylinder towards the upstream direction.

Two extreme values can be identified in all three numerical experiments in Figure 16a—a peak whose magnitude we shall refer to as the “maximum positive force” and a trough whose magnitude we shall refer to as the “maximum negative force”. In all three numerical experiments, the horizontal force exerted on the cylinder first rises to the maximum positive force, drops to the negative positive force, and then attenuates back to zero. In each of the three cases, the maximum positive force is larger than the maximum negative forces. These features are typical for a solitary wave passing a vertical cylinder, as can be seen in Yates and Wang [6] and in the validation test results (Figure 6).

Wave runup around the cylinder changes the effective area used in the calculations of wave forces. At two instants, $t=-0.1$ s, approximately when the maximum positive force occurs, and $t=0.2$ s, approximately when the maximum negative force occurs, wave runup around the cylinder is plotted in Figure 17. The polar coordinate system is used in the figure, such that $\theta =0^{\circ }$ points in the downstream direction and $\theta ={180}^{\circ }$ points in the upstream direction. It can be seen that the difference between the runup results is small. Therefore, the slight difference in the wave runup among the three numerical experiments does not explain the differences in the horizontal wave forces observed in Figure 16a.

The remaining factors causing the differences in the horizontal wave forces are then the existence of a bottom slot and the fill material placed in the slot. From Figure 16a, it can be seen that the existence of a bottom slot (i.e., test1 and test2) increases both the maximum positive force and the maximum negative force. If the slot is filled with viscous mud (i.e., test1), the maximum positive force becomes the largest of the three cases; if the slot is filled with water instead (i.e., test2), the maximum negative force becomes the largest of the three cases. A breakdown of the maximum forces is provided in

Table 3. Maximum horizontal forces exerted on the cylinder.

Name	Test1	Test2	Test3
Total maximum positive force	$1.0486$ N	$0.8613$ N	$0.8056$ N
Maximum positive force from $z=0$ to water surface	$0.7760$ N	$0.6461$ N	$0.8056$ N
Maximum positive force from $z=-0.02$ m to $z=0$ m	$0.2726$ N	$0.2152$ N	N/A
Total maximum negative force	$0.6171$ N	$0.8209$ N	$0.4949$ N
Maximum negative force from $z=0$ to water surface	$0.4274$ N	$0.5868$ N	$0.4949$ N
Maximum negative force from $z=-0.02$ m to $z=0$ m	$0.1897$ N	$0.2341$ N	N/A

.

Comparing the test2 results with the test3 results in Figure 16a reveals that both the maximum positive force and the maximum negative force are larger in test2, and the difference in the maximum negative force is more pronounced. The breakdown in Table 3 shows that the contribution from the above-slot region ($z\ge 0$ m) to the maximum positive force in test2 is actually smaller than that in test3. However, when the contribution from the slot region ($-0.02\le z<0$ m) is added, the total maximum positive force in test2 becomes larger than that in test3. In contrast, the contribution from the above-slot region to the maximum negative force in test2 is already larger than that in test3. When the contribution from the slot region is added, the total maximum negative force in test2 becomes even larger.

Comparing the test2 results with the test1 results in Figure 16a reveals that changing the fill material in the slot from water to viscous mud makes the maximum positive force larger but makes the maximum negative force smaller. The breakdown in Table 3 shows that both the contribution from the above-slot region and the contribution from the slot region to the maximum positive force in test1 are larger than those in test2; similarly both the contribution from the above-slot region and the contribution from the slot region to the maximum negative force in test1 are smaller than those in test2.

Comparing the test1 results with the test3 results in Figure 16a reveals that both the maximum positive force and the maximum negative force are larger in test1, and the difference in the maximum positive force is more pronounced. The breakdown in Table 3 shows that the contribution from the above-slot region to the maximum positive force and the maximum negative force in test1 are both smaller than those in test3. However, when the contribution from the slot region is added, both the total maximum positive force and the total maximum negative force in test1 become larger than that in test3. Therefore, although a muddy seabed is capable of attenuating wave energy, the wave energy transferred to the mud may induce mud flow, which in turn may increase the total loading on a marine structure.

To summarize, comparing the test2 results with the test3 results shows that the existence of a bottom slot increases the maximum horizontal forces exerted on the cylinder; comparing the test2 results with the test1 results shows that a bottom slot filled with mud increases the maximum positive force but decreases the maximum negative force, as compared to a bottom slot filled with water; comparing the test1 results with the test3 results shows that although a viscous mud bed attenuates wave energy, it nonetheless increases the total horizontal forces exerted on the cylinder due to the wave-induced mud flow.

Lastly, we examine the toppling moments exerted on the cylinder. The time series is shown in Figure 16b and a breakdown is provided in

Table 4. Maximum toppling moments exerted on the cylinder.

Name	Test1	Test2	Test3
Total maximum positive moment about $z=-0.02$ m	$0.0386$ N·m	$0.0323$ N·m	N/A
Maximum positive moment about $z=0$ m, calculated only from $z=0$ m to water surface	$0.0203$ N·m	$0.0172$ N·m	$0.0213$ N·m
Total maximum negative moment about $z=-0.02$ m	$0.0199$ N·m	$0.0272$ N·m	N/A
Maximum negative moment about $z=0$ m, calculated only from $z=0$ m to water surface	$0.0094$ N·m	$0.0131$ N·m	$0.0109$ N·m

. The overall trends of the toppling moments are similar to those of the horizontal wave forces—both the maximum positive moment and the maximum negative moment increase if a bottom slot is considered, and a bottom slot filled with viscous mud results in the largest maximum toppling moment. The discussions on the wave forces also apply to the toppling moments. In terms of the toppling moments, the differences between the results from the three numerical experiments become even more pronounced.

6. Concluding Remarks

This study constructed and validated a three-dimensional numerical wave flume to investigate solitary wave interaction with a vertical cylinder over a highly viscous Newtonian fluid mud bed. Numerical experiments were performed to investigate the water surface displacement, the mud surface displacement, the velocity field, and the forces and moments exerted on the vertical cylinder in the wave–cylinder–mud problem. The main findings include:

• The mud surface deformation is one order of magnitude smaller than the water surface deformation and their deformation patterns are different.
• The scour pattern caused by a solitary wave on a viscous mud bed is different from that commonly observed on a sandy seabed and warrants future investigation.
• In contrast to existing theories assuming a water–solid boundary, the existence of a mud bed instead of a solid boundary appears to have increased the thickness of the water boundary layer near the water–mud interface. Subsequently, the water flow reversal phenomenon occurring inside the boundary layer became more prominent. A more elaborate boundary layer flow analysis is needed.
• The existence of a mud bed increases the total loading and the toppling moment on the cylinder due to wave-induced mud flow.

To the best of the authors’ knowledge, no previous studies had been performed on solitary wave interaction with a vertical cylinder over a viscous mud bed, even though such a problem is of practical relevance in the design of marine structures over a muddy seabed. Therefore, this study is intended as a preliminary investigation on this subject whose main objective is to explore the underlying physical processes and to pave the way for future investigations.

Presently, scaling is based purely on existing laboratory experiments and this study is qualitative in nature. To widen the range of applicability and provide quantitative findings, more numerical experiments along with laboratory benchmark tests should be conducted to examine the sensitivity of the present findings to scaling effects, material properties of the viscous mud, the rheological model assumed for the mud, the dimensions of the mud slot, the dimensions of the cylinder, the incident wave conditions, etc. Of course, detailed information on the rheological behaviors of real mud, however difficult it may be to acquire, is essential to extending the findings to more realistic conditions. It is hoped that the findings of this study may eventually benefit the foundation design of marine structures in clay or muddy soils.