A NUMERICAL MODEL FOR WAVE PROPAGATION OVER MUDDY SLOPE

A numerical model for the interaction between waves and muddy seabed is developed, in which the motion of the movable mud and the motion of water are solved simultaneously. The governing equations for both water and the mud are the continuity equation and the equations of motion for incompressible fluids. Water is treated as a Newtonian fluid, while a visco-elastic-plastic model is used to describe the rheology of the mud. Both the interface between water and the mud and the free water surface are traced by the VOF (Volume of Fluid) method. The numerical method is based on the well-known SMAC method. The numerical model is applied to simulate wave propagation over a muddy slope, and the numerical results are in reasonable agreement with the experimental data. The present model is proved better performance than the traditional analytic model in case that topography change is not negligible.


INTRODUCTION
The existence of fluid mud has been widely reported at muddy coasts or estuaries.Fluid mud can be easily transported by gravity and by waves.It may lead to severe siltation if a navigation channel is located in such an area.As waves pass over, the oscillatory motion of the fluid mud usually causes a significant dissipation of the wave energy.Therefore, the behavior of wave propagating over a muddy bed is always one of most interesting topics in the studies of the interaction between waves and muddy bed.
During the past decades, a number of analytic studies have been carried out on the interaction between waves and a muddy bed.In the pioneering work of Gade (1958), attenuation of shallow water waves was studied based on a two-layer model, in which the upper layer, or the water, was treated as an inviscid fluid, while the lower layer, or the mud, was assumed to be a viscous fluid.Dalrymple and Liu (1978) improved the two-layer model of Gade (1958) and considered the viscosity of water.Attenuation of solitary waves and cnoidal waves were studied by Jiang and Zhao (1989) and Jiang et al. (1990).Maa and Mehta (1990) argued that the variation of the mud properties in the vertical direction can not be neglected and developed a multi-layer model.Considering the movable mud layer on seabed is thin, a two-fluid Stokes boundary layer model was developed by Ng (2000).
Study on wave attenuation due to the response of muddy seabed depends to a very large extent on the appropriate description of the rheological properties of the mud.However, the rheological property of the mud in nature has been known to be rather complicated, which depends on many factors, such as the mud density, the mineral constituents and the grain composition of the mud, the type and the concentration of ion in water, and so on.According to previous achievements, the viscous behavior of the mud is believed to be most important in wave energy dissipation.Many classical studies were based on the Newtonian viscous model (Gade 1958;Dalrymple and Liu 1978;Sakakiyama and Bijker 1989).In addition to viscosity, elasticity and plasticity also play important roles in a variety of situations.Owing to this fact, a number of more advanced rheology models were adopted to study wave attenuation due to a muddy seabed, including the visco-elastic model of the Voigt type (Macpherson 1980;Maa and Mehta 1988), the Bingham model (Mei and Liu 1987), the power-law pseudo plastic model (Huang et al. 2006) and the visco-elatic-plastic model (Shibayama and An 1993).
The studies mentioned above are all based on analytic methods.Those analytic studies on the interaction between waves and muddy bed provided a lot of useful results.Recently, by using the known analytic solution of wave dissipation over muddy bed, researchers attempted to modify the wave energy balance model to make it applicable to muddy coasts (Hall andOveisy 2007 ； Winterwerp et al. 2007;Niu and Yu 2008a).Unfortunately, all these analytic studies were limited on horizontal bed, which can not represent the real situations where topography change is not negligible.Therefore, development of a numerical method for the interaction between waves and muddy bed may be of great importance.
In this study, an effective numerical model is presented to simulate wave propagation over a region where a layer of the movable mud lies on a fixed bed, based on directly solving the continuity equation and the equations of motion for incompressible continuum.First, the governing equations and In this study, the fluid mud is considered to be a non-Newtonian fluid.The governing equations for both water and the mud are the continuity equation and the equations of motion for incompressible fluids, which means the motion of the mud and the motion of water could be described by a unified model.In this model, water is treated as a Newtonian fluid, while a visco-elastic-plastic model proposed by Niu and Yu (2008b) is used to represent the rheology of the mud.The visco-elastic-plastic model was proved of good ability in describing the rheology of the mud and modeling wave attenuation due to the response of the muddy bed.
The continuity equation and the equations of motion for incompressible continuum are where, i x are Cartesian coordinates and the indices , 1 ,2 i j = , which follow the summation convention, are used to represent the components of a vector or a tensor, t is the time, i u are the velocity components, i f are the body forces; ρ is the density, p is the pressure, and ji τ are the deviatoric stresses.To complete the formulation, a constitutive equation, which is expected to represent the physical properties of the continuum, must be added.Since water is treated as Newtonian fluid, so we have where 0 μ is the viscosity of water, which is a summation of the molecular viscosity and the eddy viscosity.The eddy viscosity is estimated by the k-ε model.And is the rate-of-strain tensor.
The rheology of the mud is complex and varies with its density and many other factors, while in this study, a visco-elastic-plastic model is chosen and the constitutive equation can be written as where y τ is the yield stress of the mud, m μ is the viscosity of the mud, G is the elastic modulus of the mud, and 2 J is the second invariant of the deviatoric stress tensor.The visco-elastic-plastic model represents a soft seabed which is essentially elastic before yield and visco-elastic after yield.

Surface and interface
The problem of interest in the present study is shown in Fig. 1.The subdomain occupied by water and that occupied by the mud can be uniformly described by Eqs. ( 1) and ( 2).In the subdomain occupied by water, the stresses are determined by Eq. (3), while in the subdomain occupied by the mud, the stresses are determined by Eq. ( 5).Therefore, the position of the movable interface as well as the free water surface should all be known at each time step to make sure the behavior of each point in the computational domain described by the correct governing equation.
In this study, both the free water surface and the interface between water and mud are described by virtue of the VOF function.In a two dimensional problem in the vertical plane, an indicator ( ) ( ) where F represents both 1 F and 2 F ; x and z are the horizontal and the vertical coordinates, u and w are the horizontal and vertical velocity components.

Boundary conditions
In order to solve the governing equations, boundary conditions for velocity and stress need to be imposed.The boundaries include the wave generating boundary, the non-reflective boundary, the rigid boundary at the fixed bed and the free surface boundary.
The left boundary of the domain acts as a wave generator so that the instantaneous water depth and the velocity there can be provided based on a wave theory by neglecting the response of the seabed, such as the theory of fifth order Stokes wave.
The right boundary is assumed to be non-reflective.In fact, the following radiation conditions are specified where η is elevation of the free water surface, and C is the phase velocity of the wave at the boundary.Since the wave is not necessarily regular at the boundary because of the complicated response of the seabed, and also because of the difficulties in estimation of the wave velocity, a subregion with artificial dissipation, called as sponge layer, is added in front of the right boundary to enhance the nonreflective boundary conditions.The artificial dissipation is realized by adding a resistance term i Du in the momentum equations, where the dissipation function D is expressed by ( ) In which, g is the gravity acceleration, h is water depth, l is length of dissipation region, x 0 is the start point of the dissipation region, n is the power of the dissipation function, θ is a non-dimensional coefficient.
The no-slip condition is applied to the fixed bed contacting with the mud, while the slip condition is applied to the fixed bed contacting with water.The no-slip condition requires that both the horizontal and vertical velocity components are equal to zero on the bottom boundary.And the slip condition requires that the velocity component normal to the fixed bed is zero while the normal derivative of the velocity component tangential to the fixed bed is zero.The pressure at the free water surface is the atmospheric pressure and equals to zero.Since the stresses and the velocity components are all continuous throughout the domain, there is no need to establish a dynamic boundary condition at the interface between water and the mud.

NUMERICAL METHOD
The motion of the muddy bed and the water are absolutely coupled.At each time step, the incident wave boundary condition is calculated first, and then the velocity components, pressure and stresses governed by Eqs.(1) to ( 5), the k-ε equations and then the VOF equations for the water surface and the water-mud interface are solved in turn, finally the computational domain and the physical parameters at each cell are renewed for the next time step.
The numerical method employed in this study for solutions of Eqs.
(1) to ( 5) is essentially SMAC, originally proposed by Amsden and Harlow (1970).Specifically, the following explicit scheme of Eq. ( 2) is used to provide a primary estimation of the velocity at ( ) is the time increment of the computational step): where the superscript n denotes the known values at t n t = Δ , i u is the primary estimation of the velocity, and To improve the values of the velocity at ( ) + Δ so that the continuity equation is exactly satisfied, we introduce a potential function Ψ and assume that Consequently, the potential function must satisfy the following Poisson equation: Taking the continuity equation into consideration, it can also be shown that Ψ can eventually be related to the pressure through ( ) The algorithm for time stepping now becomes straightforward.Based on the known values of the velocity components, the pressure, and the stresses, a primary estimation of the velocity components at the computational time step can be obtained with Eq. ( 10).The function Ψ is determined by solving the Poisson equation ( 13).The values of the velocity components at the computational time step can then be improved according to Eq. ( 12) and the pressure at the computational time step can be determined from Eq. ( 14).Once the velocity components are known, the stresses can be readily computed by using Eq. ( 3) or a simple finite difference scheme of Eq. ( 5).
The constitutive equations for both the Newtonian fluid and the visco-elastic-plastic fluid can be written in uniform expression as following.
For Newtonian fluid, 0 , and The stresses are approximated by where the physical parameters 0 1 1 , , α α β defined at the center of a cell.In the cell crossed by the water-mud interface, the physical parameters are obtained by weighted averaging.
Spatial discretization in the present study is carried out over a staggered mesh, as shown in Fig. 2. The shear stress is applied at the corner of a cell, and the normal stresses are applied at the center of a cell.The governing equations are approximated by finite differences.The advection-diffusion equations are solved with a hybrid method combining upwind scheme and second-order central scheme.The Poisson equation is discretized with the classical five-point difference scheme.
The free water surface and the moving water-mud interfaces are treated by the VOF method.The advection equation for the VOF function is solved with the donor-acceptor method so that the discontinuity of this function, from which the free water surface and the interface between water and the mud are identified, can be kept traceable.Details of the implementation can be found in Hirt and Nichols' (1981) paper.

Wave model validation
In order to validate the wave model, the model is applied to modeling wave over a submerged bar and the simulated results are compared with the observations of experiment carried out by Beji and Battjes (1993).This experiment has been frequently used in validating wave model (Beji and Battjes 1994;Gobbi and Kirby 1999).
Beji and Battjes' experiment setup is shown in Fig. 3, which is conducted in a wave flume of 37.7m in length, 0.8m in width, 0.75m in height.The still water depth is 0.4m.The top of submerged bar is 0.3m above the bed of the flume and 2.0m in length.The front slope of the submerged bar is 1:20, and its foot is located 6.0m far from the wave generator.The back slope of the submerge bar is 1:10.Wave heights are measured at 11 points shown in Fig. 3.The incident wave height 0.02m H = and wave period 2.02s T = is adopted.The computational domain is 30m long and 0.5m high, the sponge layer in front of the left boundary is 7m long.Even mesh is used, and the mesh size is 0.04m x Δ = , 0.005m z Δ = . Fig. 4 shows the time series of the water surface elevation at each measuring point.In the figures, solid line is the experimental data, and the dash line is the numerical results.As we can see, the numerical results agree well with the experimental data.Only at 5.7m x = , a large phase lag between the numerical result and the measured data is found.As Gobbi and Kirby(1999) said, this may be a mistake in measurement, so the result at 5.7m x = is not given out.
It's shown that the wave shape is gradually changing as wave is passing over the submerged bar.In the front slope, as water depth becomes shallower, wave height is increased, wave crest becomes sharper, wave trough becomes flatter and the nonlinearity of wave is enhanced.And due to the nonlinear effect, the wave energy is transferred from low frequency to high frequency.When wave pass over the top of the bar, wave nonlinear effect becomes weak as water depth becomes larger, then the high frequency wave forced in the wave group is released.Because wave phase velocity varies with frequency in deep water, the high frequency wave propagates slower and is gradually separated from the main wave as a secondary wave.The complex wave phenomenon is well simulated by the numerical model, which proves that the wave model is efficient enough for the further study.

Wave over muddy slope
In this chapter, the present model is applied to simulate wave propagation over a more complex bed, shown in Fig. 5, and compared with the observation of Soltanpour's (1999) experiment.His experiment was carried out in a wave flume at Hydraulic laboratory of Yokohama National University.The wave flume is 17.0m in length, 0.60m in width, 0.55m in height.The still water depth is 0.35m.The slope is 1:20, and its foot is set to be the origin of coordinates which is 4.7m away from the wave generator.A tilting box filled with mud, about 1.0m long and 0.2m deep, was fixed at the middle part of a slope, formed a test section from x=4.0m to 5.0m.The water content of mud is around 150%.The incident wave period is fixed to be 1.06s, and three run cases with different wave heights are carried out, whose incident wave heights are 1.73cm, 2.37cm, 7.13cm, respectively.Wave heights at four locations over the test section were recorded.The computational domain is chosen to be 13.7m in length and 0.5m in height.The rheological parameters of the mud for simulation are chosen as 180Pa G = , =10Pa s m μ ⋅ , y 7.0Pa τ = , refering to the mud rheology experiment carried out by An (1993) and Tsuruya et al. (1987).Fig. 6 shows the computed instantaneous water surface elevation and water-mud interface.As we can see, this model is efficient in simulating the wave shoaling, reflection and attenuation due to mud motion, as well as the motion of fluid mud.The mean wave height is calculated during 10 wave periods after the wave reach a steady condition, shown in Fig. 7, Fig. 8 and Fig. 9.The solid line stands for the numerical result, while the dashed line represents the result of Soltanpour (1999) by a simplified analytic model, and the dots are the measured wave height.It is evident that the present results have a better agreement with the measured data than the method of Soltanpour (1999).In the two cases that the incident wave height is smaller, the computed results at the position of wave gauge No.3 and No.4 are smaller than the measured data, which may caused by overestimation of the mud viscosity or underestimation of the yield stress of the mud.
Assuming that no mud is filled in the test section, the computed wave height is plotted by the dashdotted line in Fig. 10.Comparing with the result in case that the box is full of mud, shown by the solid line in Fig. 10, it's clear that wave height greatly decays when mud exist.
Soltanpour observed that the wave height at wave gauge No.2 was larger than that at wave gauge No.1.He ascribed the phenomena to the weakening of wave decay rate due to the fixed boundary of trench.However, we found that the effect of reflection due to the boundary of trench may even play a more important role in this phenomenon.The existence of such a reflection is clearly shown by the periodic change of the wave height in the direction of wave propagation.In Fig. 10, the dot line shows the computed wave height in case that wave propagation over a smooth slope.It can be seen that the wave reflection in case of the smooth slope bed is not notable, but in the other two cases that the trench exists, the abrupt change of the bed topography results in significant wave reflection.

CONCLUSION AND REMARKS
A numerical model for the interaction between waves and muddy bed has been developed, which ignores the exchange of sediment between the water layer and the mud layer and treat water and the movable mud as two unmixable fluids.The interface between water and the mud is tracked by the VOF method.The motion of water and the motion of the mud are obtained by numerically solving the continuity equation and the equations of motion for incompressible fluids.The numerical method directly extended from the SMAC method is adaptive to a variety of mud rheology models.
The model can be used to simulate the wave propagation over a fixed bed with a soft mud layer.To validate the wave model, the model was applied to simulate wave passing over a submerged bar, and compared with the observations of experiment carried out by Beji and Battjes (1993).It is proved efficient enough to simulate the complex wave phenomenon.
Further, the model was applied to simulate wave propagating over a slope with a section filled with mud and reasonable results were obtained.The computed wave heights were compared with the observations of Soltanpour's (1999) experiment, and it is evident that the numerical model is better than the method of Soltanpour (1999).Above all the present model is capable to deal with more complex wave phenomena, such as nonlinearity and reflection, which might be of significance for the study of the interaction between waves and muddy bed.Furthermore, this model can be freely extended to other mud rheology model and nonlinear wave condition.
t is introduced to define the free water surface. 1 1 F = implies a position occupied by fluid, while 1 0 F = implies a position occupied by air.A similar indicator t is introduced to define the interface between water and the mud. 2 1 F = implies a position occupied by the mud, while 2 0 F = implies a position occupied by water or air.Both 1 F and 2 F are governed by the advection equation as following.
Figure 1.The computational domain

Figure 2 .
Figure 2. Definition of the computational grid

Figure 4 .
Figure 4. Compare of the numerical results and measured data on Water surface elevation (measured data; ---numerical results)