Stability of mud-water interface under long surface waves

Fluid mud often exists in coastal areas with an interface separating it from its upper water layer. When a surface wave propagates over a bed covered with water and ﬂuid mud, it will cause an interfacial wave of the mud-water interface, which damps the surface wave and results in mass transport of ﬂuid mud. Most researches about wave attenuation and mass transport of ﬂuid mud are based on the assumption that the mud-water interface is unbroken. This assumption excludes the breaking interfacial waves that are known as an important mechanism responsible for mass and momentum transport between the two ﬂuids. When the surface wave is long, its velocity ﬁeld, which also serves as basic ﬂows, may be susceptible to the Kelvin-Helmholtz (K-H) instability if the shears at the interface are strong enough. In the present paper, the critical conditions for the K-H instability to occur for the mud-water interface is investigated via linear stability analysis and numerical simulation. It is found that, for a K-H instability to occur, the Stokes boundary layer thickness induced by a surface wave must be large enough to penetrate the ﬂuid mud layer and produce a strong shear at the interface. Meanwhile, a critical condition is found for a long surface wave to cause breakup of mud-water interface through K-H instability. This is practically instructive for waterway and harbor construction and protection because it predicts that a thicker mud layer is harder to be taken away by a surface wave.


INTRODUCTION
Fluid mud is typically a mixture of fine cohesive clay particles, organic matters and water, and exhibits Newtonianlike behavior. 1,2 In many real situations, the water body is stratified into two layers with a sharp interface, which is for the first time termed lutocline by Parker and Kirby. 3 The typical density of fluid mud ranges from 1050 ∼ 1200kg/m 3 and its dynamic viscosity varies from ten to thousands times that of water. When surface waves travel over a fluid mud bed, interfacial waves occur correspondingly at mud-water interface, and the energy of the surface waves dissipates due to the motion of mud. 4 The interfacial waves will cause mass transport in fluid mud layer. 5 This process is known as wave-mud interaction and has been a topic of interest 6-8 for long.
Fluid mud is often mobile in coastal areas. It is now and then suspended or deposited in channels or harbors. The bathymetric change resulted from mud transport often reduces navigability. In this case, dredging is often necessary to maintain navigation depth, which is very costly. Therefore, understanding mud transport is of significant importance for waterway construction and maintenance.
The theoretical investigation of the wave-mud interaction was first undertaken by Gade,4 who considered a two-layer fluid system where the mud layer was treated as a Newtonian fluid and the water as inviscid fluid. By adopting a shallow water assumption, Gade derived a dispersion relation, giving the wave height attenuation rate of surface waves propagating over fluid mud bed. Since the original work by Gade, the two-layer fluid system has been extended by many authors to different cases.

ARTICLE
scitation.org/journal/adv Dalrymple and Liu 8 extended Gade's two-layer fluid model by assuming the two fluids viscous and adopting the linearized Navier-Stokes equations to model small amplitude wave training with arbitrary wave length propagating in two layers of fluids. In the frame of linear waves, all variables, including the surface and interfacial velocities and pressures, were assumed to be periodical, and the surface and interfacial waves were assumed to have the same wavelengths as well as decay rates. Then, two fourth-order ordinary differential equations were obtained for the magnitude of vertical velocities. The complex wave number was to be determined from the boundary conditions, provided the water depth, mud thickness, and wave period are known. Once the complex wave number is obtained, the spatial wave height attenuation rate, the field velocity and the mass transport induced by the surface wave can be obtained. This procedure was followed by many others who considered a two-layer model in which water was assumed to be Newtonian fluid and mud was of different rheology. 7,[9][10][11][12] By introducing equivalent viscosities, mud with different rheology can also be dealt with as a Newtonian fluid.
In these previous studies, mud-water interface was assumed unbroken and no attention was paid on the stability of the interface. As a matter of fact, the interface may break due to Kelvin-Helmholtz (K-H) instability when the shear at the interface is strong enough. It is generally acknowledged that breaking interfacial waves are responsible for a significant portion of transport of mass and momentum. Once the interface breaks, the fluid mud is brought into the water column considerably. Therefore, it is practically necessary to study the critical condition for onset of the K-H instability of the interface. As we know, a longer surface wave exerts stronger shear to the mud-water interface. Hence, we investigate the critical conditions for the breaking of mud-water interface by long surface waves in the present paper.
The paper is organized as follows. After problem formulation, we simplify the velocity fields induced by long surface waves into oscillatory two-layer Poiseuille flows generated by oscillating pressure gradients. The critical condition for onset of the K-H instability in an oscillatory two-layer Poiseuille flow is then investigated using a combination of linear stability analysis and numerical simulations. Finally, some conclusions are drawn.

PROBLEM FORMULATION
The flow configuration is shown schematically in Figure 1, where the upper layer is water with depth h and the lower layer is fluid mud with thickness d. The densities and viscosities of water and fluid mud are denoted by ρ 1 , µ 1 and ρ 2 , µ 2, respectively. The coordinates x and y are along and perpendicular to the undisturbed interface, respectively, with the origin of y located at the interface. ξ, ηare the amplitudes of surface waves and interfacial disturbance waves, respectively.
For convenience, three dimensionless variables are introduced as follows. r = ρ 2 /ρ 1 is the density ratio; m = µ 2 /µ 1 is the viscosity ratio; n = d/h is the thickness ratio. The fluids are assumed to be incompressible Newtonian fluids, and they are immiscible, with a zero thickness horizontal interface between them. In both layers, the continuity and Navier-Stokes equations are employed to describe the motions of water and mud, where (u j , j ) are the velocities in x, y directions, j = 1, 2 denote the upper and lower layers, respectively, g is the gravitational acceleration.

OSCILLATING TWO-LAYER POISEUILLE FLOW
From the stability point of view, the method adopted by Dalrymple 8 is a spatial stability analysis with the basic flow being zero. All perturbation waves are attenuated, because they can't receive energy from the basic flow to develop instability. The interface is unbroken during its evolution, no matter how strong the shear induced by the surface wave is. This conclusion is based on assuming no interfacial disturbance is added to the interface. In fact, when the shear at the interface is strong enough, the wave-induced profile, which also serves as a basic flow, may become susceptible to the K-H instability which states that the interface will be unstable if the disturbance wave is short enough for a given shear. To illustrate this, the velocity field induced by a long linear surface wave is simplified here as an oscillatory two-layer Poiseuille flow, 13 which has very similar profile to the former one. The oscillatory two-layer Poiseuille flow has an exact analytical solution that can be readily obtained. This makes it convenient to analyze the stability and the motion of the interface. The difference between the long wave inducing velocity field and the oscillatory Poiseuille flow is that the former decays with time while the latter is periodic in time. The undecayed periodic flow is a little different from the flow induced by a long free surface wave, but it is suitable for flows induced by surface waves driven by winds in real environment. The surface wave-driven two-layer flow is stable if the oscillatory ARTICLE scitation.org/journal/adv two-layer Poiseuille flow is stable, because the latter flow has a stronger shear at the interface. We assume that the surface wave is long such that the vertical velocity component of the wave induced velocity filed is very small and the nonlinear term of horizontal momentum equation can be ignored. Thus, the horizontal momentum equations are linearized as where p 1 , p 2 are the hydrodynamic pressures in the upper and lower layers, respectively. Since vertical motion of the interface is small enough, the vertical momentum equations can be simplified as follows Hence, we have p 1 (x) = p 2 (x) near the interface. We consider an oscillatory flow with its far field freestream velocity expressed as U = U 0 sin(ωt), where U 0 is the velocity amplitude and ω is the angular frequency of the wave. From Euler equation, the stream-wise pressure gradient can be written as Substituting the pressure gradient The flow configuration considered here is shown schematically in Figure 2. The Stokes boundary layer thickness in the water layer induced by oscillation is δ 1 = 2ν 1 /ω, and the Stokes boundary layer thickness in the fluid mud is The velocity profile of the upper layer in two-layer Poiseuille flow is nearly uniform in the vertical direction if the thickness of upper layer d is far larger than the Stokes boundary layer thickness δ 1 near the upper boundary. Since only the short wave interfacial instability, i.e., K-H instability, is concerned, a computational thickness d 1 ≤ 20πδ 1 will be enough for theoretical analyses and numerical simulations because the influence depth of oscillation is only 2πδ 1 . The

FIG. 2.
The oscillatory two-layer flow model. basic flow for the modified oscillating two-layer Poinseuille flow can be obtained by applying a slip boundary condition at the upper layer and the no-slip boundary condition at the lower rigid wall The equations in (8) can be nondimensionalized by adopting U 0 , d 1 as the scales and written in complex form The equations have the following complex form of solutions where α 1 = (1 + i)β, α 2 = √ r/mα 1 , and A 1 , B 1 , A 2 , B 2 are integration constant determined by imposing the boundary and interfacial conditions as follows which represent slip velocity condition at the upper layer, ∂y (1, t) = 0, no-slip velocity condition at the lower rigid wall, U 2 (−n, t) = 0, continuity of velocity and shear stress at the interface y = 0, U 1 = U 2 , ∂U 1 /∂y = m∂U 2 /∂y, respectively. Substituting the integration constants obtained from the simultaneous Eqs. (11) into Eqs. (10), we have ARTICLE scitation.org/journal/adv It can be derived from Eqs. (12) that the velocity difference at the interface is: (13) The velocity profile calculated by taking the real part of Eq. (13) is shown in Figure 3, where T is the wave period, a is the wave amplitude. The maximum velocity difference at the interface during one period as a function of d 2 /δ 2 , which is proportional to n if β, m, r are fixed, is shown in Figure 4. It can be seen that the maximum velocity difference decreases with d 2 /δ 2 if d 2 /δ 2 < 3, and approaches to a constant when d 2 /δ 2 ≥ 2π. This variation tendency of the maximum velocity difference is unfavorable for a K-H instability development, which requires a sufficiently strong shear at the interface.

ONSET OF THE K-H INSTABILITY
The K-H instability is associated with steady parallel shear flow. It is caused by the shear across the interface between two fluids, and its growth rate increases with the shear  0, t))) at the interface is more unstable than a two-fluid flow with a periodic shear strength of Real(U 1 (1, t) − U 2 (0, t)) at the interface.
According to the linear theory, the dimensional critical velocity for onset of an inviscid K-H instability is: 14 where α is the dimensional wave number, and l is the disturbance wavelength. It can be inferred from Eq. (14) that the critical velocity of the K-H instability decreases with decreasing disturbance wavelength l, suggesting that a short disturbance wave is more unstable than a long one for a given velocity difference at the interface. Figure 5 shows the critical velocity as a function of the disturbance wavelength. The wave parameters and the physical properties of the fluids have been chosen such that the flow is as less stable as possible. The horizontal line in Figure 5(a) represents the maximum velocity difference between the two fluids at their interface during one period for d 2 /δ 2 = 2π. It can be inferred that, for an inviscid K-H instability to occur, the disturbance wavelength l must be less than about 0.35δ 2 . This conclusion is based on the assumption that viscosity and the boundary layer thickness near the interface have no dissipation effects on K-H instability. For a viscous K-H instability to occur, a shorter disturbance wavelength is required to overcome the viscosity dissipation without consideration of the boundary layer thickness near the interface. On the other hand, the disturbance may be attenuated due to viscosity effect if the wavelength is further decreased and the effect of the boundary layer is considered. Thus, a decrease of disturbance wavelength does not always result in a K-H instability if d 2 /δ 2 ≥ 2π. When d 2 /δ 2 becomes small, as shown in Figure 5(b), the velocity difference at the interface is quite larger than the inviscid K-H critical velocity for disturbance with wavelength equal to 0.35δ 2 . Hence, a viscous K-H instability becomes quite likely to occur. When other parameter spaces of fluid mud (ρ = 1050 ∼ 1200kg/m 3 , µ = 0.05 ∼ 10Pa · s) are considered, similar conclusions can be drawn.
To complete the analysis, we carry out a viscous analysis based on steady two-layer Orr-Sommerfeld equations governing the stability of the basic flow represented by Eq. (12) with parameters shown in Figure 5(a). Once the basic flow is determined, a systematic viscous linear stability analysis based on unsteady two-layer Orr-Sommerfeld equations with periodic coefficients is possible. According to the numerical tests carried out by Talib, 15-17 a viscous linear analyses based on the unsteady two-layer Orr-Sommerfeld equations and Floquet theory with a wide space range of parameters will cost high computational price. We can make a compromise between inviscid analyses and viscous analyses based on the unsteady Orr-Sommerfeld equations and Floquet theory. A viscous analysis with a quasi-steady basic flow is adopted instead. The quasi-periodicity approach is valid in detecting K-H instability although it cannot detect the parametric-instability that may be present in periodic problems. Fortunately, the parametric-instabilities seldom leads to breakup of the interface. [18][19][20] Thus, a quasi-steady analysis is sufficient for the present problem.
The basic flow used here is given by

ARTICLE scitation.org/journal/adv
The subscript a denotes an average of the periodic flow. The difference between the average flow and the periodic flow can be seen in Figure 6. The former doesn't change with time while the latter undergoes a complex variation during one period. The average flow rather than the periodic flow is used as the basic flow in stability analysis because the latter one is inconvenient for stability analysis and the abrupt jump of velocity near the interface is still represented by the former one, although the shear near the interface given by the latter one at some moment among one period is far stronger than that of the former one. It is noted that for d 2 = 2πδ 2 , the basic flows shown in Figure 6(a) and (b) are the same.
The viscous stability of the average basic flow is governed by the steady two-layer Orr-Sommerfeld equations and the corresponding boundary conditions as follows 21 Re ikU 1 ( where φ j (j = 1,2) represents the disturbance wave amplitude, n the wave amplitude of the interface, Fr = gd 1 /U 2 0 the Froude number, We the Weber number which is taken as zero to make the flow more unstable. The stability problem described by the governing differential equations together with the conditions at the solid boundaries and the interface represents a generalized eigenvalue problem, in which the wave speed c is the eigenvalue. A Chebyshev collocation method 22 is used to evaluate the eigenvalue problem. Once c is obtained, the growth rate ω i = c i k can be readily obtained. If ω i > 0, the flow is unstable, otherwise it is stable or neutrally stable. Figure 7 displays the growth rate as a function of disturbance wavelength at different water depths. It should be noted that δ 2 varies with water depth. It is seen from Figure 7(a) that all disturbances at different water depths decay for d 2 /δ 2 = 2π. The maximum attenuation rate occurs at a critical wavelength which is much smaller than δ 2 . A decrease of wavelength larger than the critical wavelength results in a larger attenuation rate, and an increase of water depth result in a smaller attenuation rate. On the other hand, if fluid mud thickness is decreased to d 2 /δ 2 = 1, the basic flow represented by Eq. (12) is unstable to sufficiently short wavelength disturbances for all water depths. An increase of water depth results in a higher growth rate and narrower range of unstable wavelength. The viscous analysis shown in Figure 7(b) predicts a much shorter wavelength compared with the inviscid analysis shown in Figure 5(b). This is because the basic flow used in inviscid analysis is much stronger than that of the average basic flow adopted here.
Therefore, d 2 /δ 2 plays a key role in determining the stability of the interface. When the fluid mud layer is very thick, for an inviscid K-H instability to occur, the disturbance must be very short and the viscosity may stabilize the interface. On the other hand, according to the viscous analyses based on the AIP Advances average basic flow, the K-H instability is absent at the interface, and all short wavelength disturbances decay. Since the most unstable basic flow has been adopted in calculations, it can be declared that the interface is always stable at all water depths as long as the fluid mud layer is far larger than the Stokes boundary layer thickness in the fluid mud layer. When the fluid mud layer is thin compared with the Stokes boundary thickness, both inviscid and viscous analyses predict that the interface becomes susceptible to the K-H instability. In the following section, we will carry out numerical simulations to investigate the effects of viscosity and verify the effect of fluid mud thickness on the stability of the interface at the same time.

NUMERICAL SIMULATIONS
The flow configuration is shown schematically in Figure 2. Eqs. (1-3) are discretized using uniform control volumes and staggered variable arrangement is employed. The PISO algorithm 23 is adopted to solve momentum and continuity equations. In additional to Eqs. (1-3), the volume of fluid method(VOF) is adopted to capture the interface between different fluids. The equation for VOF is where F is the volume fraction function which represents the volume ratio occupied by the fluid mud in a cell, ρ and µ are the density and the dynamic viscosity of fluid in each cell, respectively, which are determined by the following expressions. At the bottom and top boundaries, no-slip and slip boundary conditions are imposed, respectively. Periodic boundary condition is applied in the stream-wise direction. The computational domain L x × L y is divided into M × N grids, and grid coordinates are x i = iL x /M, x j = jL y /N (where i = 0, 1 · · · M; j = 0, 1 · · · N). Grid independence tests are performed by increasing M and N until there is no qualitative difference in the flow patterns. M × N = 220 × 300 is found to be sufficient for all the cases considered below. Our code is validated with exact laminar solutions of oscillatory flow velocity profile by setting r = 1, m = 1 in Eq. (12). Excellent agreements are obtained ( Figure 8).
Traditionally, the simulations are seeded with the eigenmode obtained from viscous stability analyses. Here, the disturbances are introduced in another manner, which is convenient to carry out analyses in both numerical simulations and lab experiments. At the initial time, the fluids are motionless, and then a disturbance wave given by a standing wave theory is introduced to the fluids to obtain an approximately initial velocity field. The standing wave theory gives 25 where η is the amplitude of the interface, and A 0 = l 20 is the corresponding maximum wave amplitude. Setting t to be zero gives an initial flow field with η = 0. Subsequently a periodic pressure gradient is applied in the stream-wise direction to drive the fluids. The disturbance will decay because of viscous dissipation if the periodic pressure gradient is absent, and will grow if the shear due to the pressure gradient is strong enough to trigger an instability. The stabilizing effect  Figure 9 shows the evolution of the disturbance with wavelength l=8d 2 . The computational box contains two consecutive waves. It is clear that the disturbance decays as time elapses and the flow is stable. This is true even if the disturbance wavelength is further decreased to l = 4δ 2 and l = 2δ 2 , as shown in Figure 10, Figure 11. Furthermore, as the disturbance wavelength gets smaller and smaller, the time for the disturbance to decay becomes shorter and shorter. This indicates that viscosity has a dissipation effect on the disturbances because shorter disturbance is less stable if the viscosity is unconsidered, as demonstrated by Eq. (14). This is consistent with the finding obtained from the above viscous stability analysis, which predicts that the decrease of disturbance wavelength larger than the critical wavelength (only one fifth or less of δ 2 ) will result in a more and more stable flow.
Hence, it can be anticipated that a further decrease of the disturbance wavelength will result in a more stable flow region and the K-H instability will be absent at the interface. Now, let us examine the effect of mud layer thickness on the stability of the interface. Figure 12 shows the snapshots of the interface when it is close to the breaking point for different mud layer thicknesses. It is clear that the thinner the mud layer is, the earlier the break begins, indicating the flow is more and more unstable. This is consistent with the finding obtained from the linear analyses. Figure 13 shows the evolution of the interface at different time with the mud layer thickness d 2 = δ 2 . It can be seen that the interface breaks due to a shear instability if the fluid mud is thin enough. One may note that the unstable wavelength exceeds the scope of unstable wavelength predicted by viscous analysis as shown in Figure 7(b). This is because the shear here is stronger than that in viscous analysis shown in Figure 7(b). It should be noted that the surface tension has been ignored and the oscillatory two-layer Poiseuille flow has been adopted to make the flow more unstable. This makes it favorable for formation of the K-H instability, but it doesn't change the fact that the fluid mud thickness plays a center role in determining the stability of the interface between fluid mud and water. Since the effect of surface tension is stabilizing, it is reasonable to anticipate that for the interface to break, the surface wave must be very long, or equivalently, the fluid mud layer must be very thin. On the other hand, the interface will be stable if the fluid mud is very thick.

CONCLUSIONS
Mud-water interface, also called lutocline that separates water and fluid mud, often exists in coastal areas. Its stability gives significant impact on the fluid mud suspension, which influences navigation depth greatly. Wave theory tells us that a longer surface wave produces stronger shear at mud-water interface. Therefore, mud-water interface can be susceptible to the K-H instability if a surface wave is long enough. We have investigated this issue in the present paper by simplifying the wave inducing velocity fields to the oscillating two-layer Poiseuille flows.
The critical conditions for the K-H instability to occur have been examined for different physical parameters of the two fluids and surface waves. The linear stability analyses show that the ratio of the mud layer thickness and the Stokes boundary layer thickness in mud layer induced by surface waves plays a central role in determining the stability of the interface. The flow will be stable if the Stokes boundary layer thickness is far smaller than the mud layer thickness, and it will become susceptible to the K-H instability if the mud layer is sufficiently thin. Numerical simulations are then performed as complementary evidences to examine the effects of disturbance wavelength and mud layer thickness on the ARTICLE scitation.org/journal/adv stability of the interface. The simulations show that a shorter disturbance wave results in a more stable region if the mud layer is far larger than the Stokes boundary layer thickness in the fluid mud layer. On the other hand, the flow becomes more unstable if the mud layer thickness decreases. The mudwater interface does break if the Stokes boundary layer thickness is large enough to penetrate the mud layer. This is consistent with the finding from our linear stability analyses. A sufficiently long surface wave, or equivalently a sufficiently thick Stokes boundary layer, will cause instability of the interface because of strong shear at the interface. This is practically instructive for waterway and harbor construction and protection.