A hybrid numerical model is introduced for simulation of cohesive sediments concentration profiles in a surf zone. For this purpose, wave height reduction must be considered, due to muddy beds and wave breaking. Models, such as Sanford and Maa’s erosion model, Krone’s sedimentation model, Tajima’s wave breaking model and the viscoelastic–plastic rheological model, are used to investigate the interaction of wave and bed and to predict the concentration profile. A splitting algorithm has been used to split the threedimensional advection–diffusion equation into a horizontal, twodimensional equation, and a vertical, onedimensional equation, due to different length scales. The onedimensional equation is discretized over a non uniform grid, and, then solved implicitly using the QUICKEST scheme (third order in time and space). The twodimensional equation is divided into two parts (advection and diffusion) and each part is separately solved at different time steps. A suitable mesh, regarding space and time intervals, is chosen for considering the stability of the present model. The computational domain extends from the shoreline to the deepwater zone. Finally, the results are analyzed and compared with experimental and field data and other models. Good agreement has been obtained with the data and other numerical models.
Hybrid model ; Concentration ; Cohesive sediment ; Split method ; Surf zone
Numerical simulation models have great potential for use in estimating the dredging effects in ports, and the maintenance of harbors and navigation channels in estuaries, etc. Cohesive sediment can also serve as a source or sink for pollutants introduced into the water column. Numerical modeling of concentration profiles for cohesive sediment in coastal zones should be carried out by considering the specific characteristics of the cohesive sediment, the interaction between wave and muddy bed, erosion and sedimentation. Muddy bed reaction on waves compared to sandy beds is quite different because, after occurrence of fluid mud in the bed, the energy depreciates and wave height is reduced. It was commonly proved that the friction coefficient of waveinduced bed shear stress for muddy bottoms is greater than for sandy bottoms [1] .
In other words, soft mud interacts with waves, resulting in attenuation of wave height due to bottom friction, percolation losses and viscous damping within the sediment. These interactive modes are also manifested as changes in wave length, water particle motion, and the elevation of the interface between the fluid and bottom sediment [1] . Wave breaking in near shore zones is the dominant phenomenon to cause long shore and on–off shore currents, which play a significant role in sediment transportation. Thus, prediction of the height and location of waves breaking on the shore is a very important issue.
The main objective of this paper is to predict the concentration profiles of cohesive sediment in a surf zone, considering Sanford and Maa’s erosion model [2] , Krone’s sedimentation model [3] , Tajima’s wave breaking model [4] and the viscoelastic–plastic rheological model [5] , [6] and [7] for muddy beds. A splitting algorithm has been used to split the threedimensional transport equation into a horizontal, twodimensional equation and a vertical, onedimensional equation, due to different length scales. The horizontal, twodimensional equation is divided into two parts, including advection and diffusion, and each part is solved separately at different time steps [8] .
The flowchart of the present model is shown in Figure 1 . By defining input parameters, the wave–mud interaction model is used to calculate the wave attenuation rate coefficient. Wave height is calculated using Tajima’s wave model, and the breaking point is determined. Then, shear stresses are computed using Tajima’s wave model. These shear stresses are compared with critical erosion and sedimentation shear stress. Consequently, the amount of erosion and sedimentation are computed. The erosion and sedimentation rates are used in the advection–diffusion equation as sink and source terms. Suspended sediment concentration is calculated by solving the advection–diffusion equation using the standard split method. Finally, the results of modeling are analyzed and compared with experimental and field data and other models.

Figure 1. Flowchart of the present model.

Rheology, defined as the science of fluid mud deformation, is an important field in the investigation of fluid mud behavior response to wave actions. The rheological property of the mud has been known to be rather complicated, which depends on many factors, such as mud density, mineral constituents, grain composition of mud, and type and concentration of ions in the water, etc. The viscoelastic–plastic model is a combination of linear viscoelastic and non linear plastic models, and is used for mud behavior modeling. In this model, fluid mud is assumed to be viscoelastic when stress is less than yield stress. The latter is considered to be viscoplastic when applied stress is greater than yield stress. The constitutive equation is expressed as [5] and [6] :

( 1) 
where and take the values 1 and 2, which correspond to and axis, respectively, is the angular frequency of the wave, is the apparent viscosity, is the deviator part of the stress tensor, is the deviator part of the strain rate tensor, is the elastic modulus, is the viscosity of mud in the viscoelastic state, is the viscosity of mud in the viscoplastic state and is the yield stress. is the objective of the deformationrate tensor and is expressed as [5] , [6] and [7] :

( 3) 
where and are the horizontal and vertical coordinates, and are velocity components in and directions, respectively. The rheological parameters are evaluated from laboratory experiments [5] and [6] . The fluid system is divided into layers, in which the water layer is represented by [9] . For a small disturbance to this system, by neglecting convective accelerations, the momentum and continuity equations can be derived from the equation of motion for a Newtonian, laminar and incompressible fluid under wave action, as [6] and [10] :

( 4) 
where the subscripts, , indicate the layers, and parameters , , and represent the time, density, kinematic viscosity of mud and dynamic pressure, respectively. The separated periodic solutions for , and are assumed as [6] and [10] :

( 7) 
where is angular frequency and is the unknown complex wave number, which can be expressed as:

( 10) 
Displacements of water surface and interfaces, , are represented as below:

( 11) 
where is the amplitude of displacement of the th layer and the water surface is expressed as . Therefore, the real part of the wave number, , gives the wave length and its imaginary part, , is the wave attenuation rate, assuming the exponential wave height decay (Section 2.3 ). Finally, the fourth order differential equation for is solved, with respect to the prevailing boundary conditions, and the solutions are obtained. The unknown constants and variables are determined from the boundary conditions at the water surface, the interfaces and the rigid bottom [6] and [10] . Finally, wave attenuation rate coefficient, , is calculated by a wave–mud interaction model.
Tajima (2004) [4] defined the equivalent linear wave with an energy flux equal to that of the actual nonlinear wave. The model uses linear wave theory to propagate waves until the waves reach breakpoint. After breaking, a broken wave dissipation term is introduced in the energy flux balance equation. The breaking criterion for a linear wave is defined as [4] :

( 12) 
where is bed slope, is the wave number at breaking point, is the wave height at breaking point, is deepwater wave length and is the wave depth at breaking point. Tajima’s formulation for broken wave energy dissipation is:

( 13) 
where is the group velocity, is the horizontal direction, is the water depth, is the wave energy, based on the recovery wave height, and is a dissipation factor. According to Tajima, is dependent on the bed slope as:

( 14) 
where , and is the slope of the mean water depth. Considering Tajima’s experimental data, which are detailed in Tajima (2004) [4] , the values of and are obtained. Surface roller energy, is calculated by:

( 15) 
Also, energy dissipation rate is calculated by:

( 16) 
where is roller crosssectional area, is wave phase velocity and is the wavelength. The volume flux due to the surface roller, which causes an increase in the return current, is:

( 17) 
Tajima assumed the following equation to calculate the energy balance for the surface roller:

( 18) 
where is the fraction of the energy of broken wave energy that goes into the surface roller (Tajima assumes ), and is proportionality constant. Tajima takes , and the complete energy balance equation for the roller is expressed as [4] :

( 19) 
The near shore mean current velocity below the trough level is given by:

( 20) 
In the above equation, the current shear stress is assumed to vary linearly in depth. Here, is the turbulent eddy viscosity, is the mean current shear stress vector, and are the bottom and trough shear stress vectors, is the trough water depth, and is the elevation above the bottom. The surface roller acts as a source of wave energy dissipated by wave breaking. Only part of this energy is immediately transferred to turbulent energy, while the rest is temporarily stored in the surface roller. The surface roller introduces a new forcing term that accounts for the observed increase in return flow velocity. Close to the bottom, shear stress is considered constant, and the turbulent eddy viscosity is assumed to vary linearly with depth. Shear stress increases because of the effect of broken waves at the points away from the bottom. In this region, the rate of increase in turbulent eddy viscosity must be faster than the previous region. The turbulent eddy viscosity is defined as:

( 21) 
where is Von Karman’s constant, is the shear stress velocity at the bottom, is the shear stress velocity at the surface, is the depth where the two eddy viscosity profiles match in Eq. (21) , and is the wave bottom boundary layer thickness, determined from Madsen’s model [11] .
The mean current velocity at the outer edge of the boundary layer is:

( 23) 
where , is the equivalent Nikuradse roughness of the bottom. The maximum combined wavecurrent bottom shear velocity is defined by [4] :

( 24) 
With:

( 25) 
where is the angle between waves and currents. The boundary layer thickness, , is given by:

( 27) 
where:

( 28) 
and is the amplitude of the nearbottom wave orbital velocity, , whose value is provided by Tajima’s model. According to Madsen (1994) [11] , wave friction factor, , can be approximate as a function of the dimensionless parameter, as [11] :

( 29) 
For the timeaverage over a wave period, the components of the wave volume flux above the trough level are:

( 30) 
where are the wave orbital velocity components in and directions, respectively, is the wave energy density, and is the angle of wave incidence. Wave radiation stresses, and , are given by:

( 31) 
The volume fluxes in and directions due to the surface roller are:

( 33) 
where and are the surface area and energy of the roller, respectively. Averaged momentum fluxes due to the surface roller, and , are expressed as [12] :

( 34) 
According to the near shore current model, surface velocity is:

( 36) 
where:

( 37) 
where is the maximum combined wavecurrent bottom shear velocity, , is the equivalent Nikuradse roughness of the bottom, and is the wave amplitude [4] . Finally, the system of equations yields:

( 42) 
Analogously, and are calculated [4] :

( 44) 
Briefly, the main variables that define the model are: wave setup, , total flux due to currents, , surface velocities, , bottom shear stresses, , and trough shear stresses, . Since there are five independent variables, it requires five boundary conditions at each of two boundaries: five at the offshore boundary and five at the onshore boundary [4] .
The attenuation of wave height on a horizontal bed is usually approximated by an exponential function [7] :

( 46) 
where is the incident wave height at , and is the wave attenuation coefficient. The wave attenuation coefficient depends on water depth, wave characteristics, the thickness of fluid mud layer and its rheological behavior. Applying the exponential wave height decay over a muddy and gentle bottom slope, the energy dissipation rate is expressed as:

( 47) 
After simplifying, the above equation is:

( 48) 
where is the wave energy per unit surface area, is the water density and is the group velocity.
Erosion is one of the major factors in sediment resuspension, sediment transport and beach deformation of cohesive shorelines. Erosion of cohesive sediments occurs when hydrodynamic erosive forces exceed gravitational, cohesive or frictional forces (surface erosion). The second condition occurs when flowinduced shear stress exceeds bed bulk shear stress, which is called mass erosion. The simple form of unlimited erosion linearly correlates as:

( 49) 
Sanford and Maa (2001) chose the linear erosion rate formula as [2] :

( 50) 
where is erosion rate, which is a function of depth and time , is erosion constant, is the bed shear stress and is the critical bed shear stress of sediment erosion. The basic assumptions were as follows: (a) Erosion constant is proportional to sediment concentration at the water–mud interface; (b) Critical shear stress increases locally at a constant rate with depth. Therefore, erosion rate depends on both time and depth [13] :

( 51) 
where is the sediment density, is the volume fraction of sediment, is a local constant, is the erosion velocity and is the rate of increase in resistance against erosion, with respect to depth, that is assumed to be locally constant. By assuming the homogeneous solution ( ), the following solution is obtained:

( 54) 
where is the critical shear stress when is first applied at . Mass erosion occurs with rapid increase in the fluid shear stress in the surf zone and where the wave breaks. The mass erosion rate per unit area may be written as [13] :

( 55) 
In the above equation, is the dry density and is the thickness of the erodible layer. In this study, the values of critical shear stress for surface and mass erosion are considered 0.5 and 1.5 Pa [2] and [13] .
Cohesive sediments tend to agglomerate together. This process is called flocculation and the resulting particle is called floc. Floc size is different, which complicates prediction of settling velocity. Settling of mud floc is one of the most important parameters in evaluating the concentration profile in marine environments. Deposition rate has the following formula [3] :

( 56) 
where is the bed shear stress, is the critical shear stress for the deposition, is the nearbed particle concentration and is the settling velocity. The near bottom concentration is dependent on depthaveraged sediment concentration, as described by Teeter (1986) [14] . In this study, the value of critical shear stress for deposition is considered as 0.2 Pa [3] and [9] .
The settling velocity of cohesive sediment depends on concentration and is usually divided into three regions [9] :
(a) Free settling

( 57) 
where is the specific gravity of sediment, is representative floc diameter and is the fluid kinematic viscosity.
(b) Flocculation settling

( 58) 
where is an empirical coefficient.
(c) Hindered settling

( 59) 
where is the value of at concentration , is the inverse of and is the concentration at . The Winterwerp’s equation can also be used [15] :

( 60) 
where is the dry density of sediments, is the total concentration and is the gelling concentration of suspended sediments.
The partial differential equation describing threedimensional advection–diffusion for suspended sediments can be written as [8] :

( 64) 
where is suspended sediments concentration, , and are velocity components in the , and directions respectively; is settling velocity of suspended sediments; is the source or sink term; and , and are diffusion coefficients in , and directions, respectively. These coefficients are equal to eddy viscosities [16] . In coastal waters, the horizontal length scale is generally much larger than the vertical scale. Therefore, Eq. (64) will split into a twodimensional horizontal equation (Eq. (65) ) that is solved horizontally and, then, the onedimensional vertical equation (Eq. (66) ) that is solved vertically [17] .

( 65) 
Also, the following boundary conditions for sediment are used [17] :

( 67) 

( 68) 

( 69) 

( 70) 
where is bed shear stress, is critical shear stress for deposition, is critical shear stress for erosion, and and are erosion and deposition rates (Sections 2.4 and 2.5 )
The horizontal twodimensional equation (65) can be converted in a pure advection system by neglecting the diffusion part (Eq. (71) ) and in a pure diffusion system by neglecting the advection part (Eq. (72) ), and can be written as [8] :

( 71) 
By considering the layerintegrated hydrodynamic model, Eq. (65) was integrated over the layers to give the following twodimensional layerintegrated equation:

( 73) 
where is depth averaged suspended sediment concentration, and are the depth averaged flow velocities in the and directions, respectively, is water layer thickness, and and are erosion (source) and deposition (sink) terms in upward and downward directions, respectively. To solve Eq. (73) , a fractional step approach, also known as standard the split method, is used [18] . In this approach, both advection and diffusion parts of the advection–diffusion equation are solved separately at each time step. The splitting approach is an accurate numerical method to solve advection and diffusion separately. A high resolution conservation algorithm for the advection part in incompressible flow was developed by Leveque (1996), who used a basic upwind method and proposed several correction terms to achieve better accuracy and stability [19] . Also, to solve the diffusion part, a semi implicit finite difference scheme is employed, such that, it can easily be converted to a completely explicit or implicit scheme. Discretization of the concentration and velocity field in the solution domain is shown in Figure 2 [8] .

Figure 2. Discretization of concentration and velocity field in solution domain [8] .

A conservative form of advection of scalar concentration or density function can be written, in general, as [8] :

( 74) 
By assuming incompressible flow, a twodimensional advection equation can be written as:

( 75) 
In discrete form and for every cell with an equal size of mesh in two directions, the following condition should satisfy:

( 77) 
where is understood that all data is at time . To solve this conservative form of the advection equation, the upwind method is used, which is based on flux calculation of the concentration at the cell interfaces. By assuming time steps, , it can be written as:

( 78) 
where represents the flux at the left interface of cell , and represents the flux at the right interface of cell . Similarly, represents the flux at the bottom interface of cell , and represents the flux at the top interface of cell . These fluxes at the cell interfaces can be calculated as:

( 79) 
Therefore, Eq. (78) can be rewritten as:

( 80) 
The modifications can be incorporated in the flux calculation of and as follows:

( 81) 
To achieve a second order accurate algorithm, a second order Lax–Wendroff method is combined with the upgrade upwind method. The Lax–Wendroff method can be expressed as Eq. (82) . It can also be rearranged by considering the upwind method combined with a correction term, as Eq. (83) :

( 82) 
Similarly, this method can be used at the interface between and by adding the following term to the interface flux. To avoid oscillation, a flux limiting factor is also introduced as [8] :

( 84) 
And, similarly, for the flux at the interface between and :

( 85) 
The flux limiting term, , is defined as:

( 86) 
Some standard limiters are minmod, superbee, van Leer and monotonized centered, as Eqs. , , and :

( 87) 
The superbee limiter provides the best results. The transverse motion of the correction wave at the interface between cells modifies the fluxes and can be calculated by the following expressions [8] :

( 91) 
Transverse motion between cells, and , modifies the fluxes, and , and can be calculated using the following expressions:

( 92) 
Also, a finite difference representation for the diffusion part can be written as follows:

( 93) 
The above equation is solved for time steps, , using an explicit finite difference scheme. Introducing a new variable, , Eq. (93) can be rewritten as:

( 94) 
where:

( 95) 
Therefore, Eq. (94) is written as:

( 96) 
where and can be calculated as [8] :

( 97) 
The onedimensional equation (66) was discretized over non uniform grid spacing in the vertical direction. It was solved implicitly using the QUICKEST (Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms) scheme and the backward or implicit Euler method. The coefficients at the resultant discretized equation are computed once and stored in arrays, and then they are used to update each constituent. The coefficients matrix has three diagonals with nonzero elements, and the values of concentrations are computed by inverting the resulting matrix and by using the Thomas tridiagonal solver. The resultant discretization equation is [17] :

( 99) 
where superscript is an index for time and the subscript denotes layer number at depth (at free surface, ). Also, other parameters are defined as:

( 100) 
where the symbol is used to denote the greater value between and , superscripts and denote the control volume face (top and bottom), is the grid Peclet number and is the mass flow rate. The Peclet number is defined as the ratio of the rate of advection of a physical quantity by the flow to the rate of diffusion of the same quantity driven by an appropriate gradient. In the context of species or mass dispersion, the Peclet number is the product of the Reynolds number and the Schmidt number.
To model wave damping, the East Bay situated in Louisiana, USA was considered. The deposited sediments from the Mississippi River are very soft, with shear strength ranges from 1.57 kPa near the watersediment interface to 2.36 kPa three meters into the sediment. The wave period was 7.75 s and the measured wave heights were 0.285 and 0.54 m at two stations, A and B (Figure 3 ). The wave height at B is set to be 0.54 m and the wave height transformations are calculated based on Tajima’s wave model. The bottom contours are shown in Figure 3 [6] .

Figure 3. Bathymetry and wave measurement stations in the East Bay, Louisiana [6] .

Soltanpour’s experimental data are used for validation of the present model. The experiments were carried out in a wave flume at the Hydraulic Model Laboratory, Department of Civil Engineering, K. N. Toosi University of Technology. These laboratory experiments were performed in a wave flume; 17 m long, 0.6 m wide and 0.55 m deep. Water depth was kept constant at 35 cm in all experiments, wave period was 1.06 s, and the wave heights were varied during the tests. The mud was made using a mixture of kaolinite and water set on the bed [5] . A concentration profile of cohesive sediments was modeled using two sets of field experiments along the Louisiana coast (US). One of these measurements was performed by Kemp (1986) [20] at 30 km of coastline, southwest of Vermilion Bay, on the eastern Louisiana Chenier Plain. The sedimentation in the area is due to the Atchafalaya River. The bed material consists of about 50%–80% clay, 20%–50% silt and zero to 14% sand. Mud depth was found 0.18 m and mud density was taken 1270 kg/m^{3} . Measured water depth varied between 0.85 and 1.05 m and water density was considered 1030 kg/m^{3} in the region. The second data was collected using WAVCIS (WaveCurrentSurge Information System) in a 20 km length of south of Marsh Island along the Louisiana coast. These data are obtained by Sheremet and Stone (2003), Sheremet et al. (2005) and Sheremet (2006), and were reported by Sorourian (2007) [13] . The study site is located along the 5m isobaths, on a sloping shelf (<0.001). Bed material is inorganic fine sediment with a median grain size of 6.34 μm. The mean significant wave height was found to be 0.177 m and the mean significant wave period was 3.17 s [13] . It is worthwhile to mention Jazayeri Shooshtari and Soltanpour’s model and experimental data, whose concentration profiles are used for validation. Experimental data were obtained by Van Rijn and Louisse (1987) and conducted on a kaolinite bed with under wave action. These experimental data were reported by Jazayeri Shooshtari and Soltanpour (2008) [21] .
Regarding the stability of the hybrid numerical model, space and time intervals must be chosen correctly. The bathymetry of the computational domain extends from shoreline to the deepwater zone. A suitable mesh was constructed for covering the study area. The spatial increments ( , for concentration profile, time step=1 s) in the wave breaking zone and ( , for concentration profile, time step=1 s) in the outer region were used to create a uniform mesh. In order to find the proper mesh size for discretization, a sensitive analysis was carried out for 600 m in a cross shore direction. The parameters, such as bed slope, deep water wave height, wave period, deep water wave angle, mud bed thickness and water content of mud, as 1%, 1.05 m, 7 s, 45°, 25 cm and 120%, respectively, were considered for sensitivity analysis. The Tajima’s wave model was run for 5 different mesh sizes. The wave heights at depth equal to 2.75 m were compared with the wave heights calculated from the exact solution (1.201 m), which is based on linear wave theory. The results are shown in Table 1 . The region in the surf zone is within 100–200 m of the shoreline, and out of this region is within 200–600 m of the shoreline.
Case  Mesh size/surf zone (m)  Mesh size/out of surf zone (m)  Wave height at point with depth = 2.75 m  Relative error (%)  Wave height at breaking point (m) 

1  1  10  1.203  0.17  1.403 
2  5  15  1.257  4.66  1.464 
3  10  20  1.291  7.49  1.502 
4  15  25  1.397  16.32  1.624 
5  20  30  1.504  25.23  1.772 
Figure 4 and Figure 5 show validation of present model by field data from Louisiana East Bay, Soltanpour’s model and its experimental data [5] . The results of comparison show that there is good agreement. Figure 4 shows wave height distribution in Louisiana East Bay with different mud bed thicknesses for the present model, Soltanpour’s model and field data [5] . As can be seen, the present model is matched well with Soltanpour’s model. For mud bed thickness equal to 50 cm, the difference in distance of the breaker line is about 5 m closer to the shoreline, and the wave height at breaking point is reduced about 4 cm relative to Soltanpour’s model. The results from field data and mentioned models show good agreement with mud bed thickness equal to 70 cm. Figure 5 shows wave height distribution in the experimental flume for the present model, Soltanpour’s model and experimental data (second case) [5] . As can be seen, the present model is matched well with Soltanpour’s model. The result of the present model shows that the difference of breaker line in distance is about 4 cm closer to the shoreline, and wave height at breaking point is about 4 mm less than Soltanpour’s model.

Figure 4. Wave height distribution in Louisiana East bay with different mud bed thicknesses (present model, Soltanpour’s model and field data [5] ).


Figure 5. Wave height distribution in experimental flume (present model, Soltanpour’s model and experimental data, second case [5] ).

Another part of the validation for the suspended sediment concentration profile was to compare the results of the present model with Sorourian’s model and field data. The field data from the Louisiana coast are comprised of one set of data by Kemp (1986) [20] and the second set of data by Sheremet (2006). Sorourian made the following assumptions in his model: (a) Linear wave theory; (b) Onedimensional and waveaveraged; and (c) Horizontal initial bed. He incorporated several theoretical concepts to enhance the accuracy of the suspended sediment concentration predictions, such as considering the electrochemical force between cohesive sediment particles, relating bed resistance to sediment characteristics, predicting bed shear stress by relating the wave friction factor to bed roughness and median grain diameter, and considering timedependent erosion rate [13] . Figure 6 and Figure 7 show the samples of comparison between concentration profiles of the present model with Sorourian’s model and field data. The results in these figures show that there is very good agreement with field data. The maximum suspended sediment concentration difference between the present model and field data is 3 g/l at depth 0.12 m (Figure 6 ) and 0.03 g/l at depth 1 m (Figure 7 ). As can be seen from these figures, the maximum difference in concentration with Sorourian’s model is 2.5 g/l at depth 0.06 m (Figure 6 ) and 0.006 g/l at depth 2 m (Figure 7 ). Due to assumptions made in Sorourian’s model, the present model predicts suspended sediment concentration more accurately.

Figure 6. Comparison of concentration profile between present model, Sorourian’s model (2007) [13] and kemp field data (1986) [20] for wave height equal to 0.117 m and wave period equal to 6.6 s after 90 min (total time=180 min, time step=1 s, mud thickness=18 cm).


Figure 7. Comparison of concentration profile between present model, Sorourian’s model (2007) [13] and Sheremet field data (2006) [13] for wave height equal to 0.177 m, wave period equal to 3.17 s and mud depth equal to 0.6 m (time: 12:59, 11/7/2003).

The results of the present model also are compared to Jazayeri Shooshtari and Soltanpour’s model and experimental data [21] , which is shown in Figure 8 . Their model is a crossshore numerical model that simulates a cohesive sediment profile under wave action. The partial differential equation of twodimensional vertical advection–diffusion is solved. A dynamic mesh is employed, where the physical space in Cartesian coordinates is replaced by a computational coordinate system. Both water level and bed surface are time dependent and the finite volume method is applied to solve the discredited equation [21] . The results of the present model show that there is very good agreement with Jazayeri Shooshtari and Soltanpour’s model and experimental data at the zone near the bed. At the upper zone, the maximum suspended sediment concentration difference between the present model and Jazayeri Shooshtari and Soltanpour’s model is 16 mg/l at depth 0.15 m. As mentioned earlier, the splitting approach is a very accurate numerical method can be used to solve advection and diffusion parts, separately [8] .

Figure 8. Comparison of concentration profile between present model, Jazayeri Shooshtari and soltanpour’s model and experimental data [21] for wave height equal to 0.03 m and wave period equal to 1.5 s after 5 min (erosion threshold shear stress equal to 0.07 Pa and deposition threshold shear stress equal to 0.04 Pa).

Comparison of the results of the present model with Sorourian’s model (2007) [13] , Jazayeri Shooshtari and Soltanpour’s model (2008) [21] , and their reported experimental data and Louisiana coast field data, confirms the accuracy and validity of the present model. Based on the present results, the advection–diffusion equation is solved by the standard split method to obtain more accurate results in comparison to experimental and field data.
By assuming mud bed thickness equal to 80 cm, the present model has been run and concentration profiles are predicted at approximate depths equal to 5 and 16 m (out of the surf zone) after 60 min of entering the wave in Louisiana East Bay [21] . The results are compared with Jazayeri Shooshtari and Soltanpour’s model results. As can be seen, the present model matches well with their model. Figure 9 (a) and (b) show the results of Jazayeri Shooshtari and Soltanpour’s model and the present model in predicting concentration profiles at these depths. To show the ability of the present model in prediction of the concentration profile within the surf zone, the model has been run for a depth equal to 1 m and mud thickness equal to 50 cm after 60 min. The result is shown in Figure 10 .

Figure 9. Comparison of concentration profile between present model with Jazayeri Shooshtari and Soltanpour’s model [21] at approximate depth equal to (a) 5 m and (b) 16 m (out of surf zone), East Bay, Louisiana.


Figure 10. Concentration profile at depth equal to 1 m (within surf zone), East Bay, Louisiana.

To introduce the application of the present model to calculating the wave attenuation coefficient, wave height, wave height to water depth ratio, bed shear stress and erosion rate, the model has been run by parameters, such as deep water wave height, wave period, deep water wave angle, water content of mud, and bed slope, as 1.05 m, 7 s, 45°, 120% and 0.01, respectively after 30 min. Variation of the wave attenuation coefficient as a function of water depth is shown in Figure 11 for a mud bed thickness equal to 70 cm. As can be observed, wave attenuation rate decreases as water depth increases. In other words, by increasing water depth, muddy bed effects are less, because deep water waves do not interact with the sea bed [1] .

Figure 11. Variation of wave attenuation coefficient as a function of water depth.

Variation of the wave attenuation coefficient as a function of wave height is shown in Figure 12 for a mud bed thickness equal to 70 cm, water depth equal to 1.2 m and different deepwater wave heights. As can be observed, wave attenuation rate increases as wave height increases. In other words, by increasing wave height, muddy bed effects are significant, because water waves have more interaction with the sea bed.

Figure 12. Variation of wave attenuation coefficient as a function of wave height.

Figure 13 show the variation of (a) wave height to water depth ratio and (b) wave height versus distance to shoreline. Based on these figures, by increasing mud bed thickness, the wave breaker line is moved to the shoreline; on the other hand, if the mud bed is too thick, the water waves may not break, due to bottom friction, percolation losses and viscous damping within the sediment.

Figure 13. Variation of (a) wave height to water depth ratio (b) wave height for different values of mud bed thickness.

Bed shear stresses were also computed over the whole domain from the wave fields through Tajima’s wave model. Figure 14 shows the variation of bed shear stress versus distance to shoreline. Maximum bed shear stress occurs at a distance of about 168 m from the shoreline (breaker line is about 172 m from shoreline). As expected, maximum bed shear stress occurs near the breaker line.

Figure 14. Variation of bed shear stress versus distance to shoreline.

Figure 15 shows variation of the surface and mass erosion rate versus distance to shoreline. As mentioned formerly, the value of critical shear stress for mass erosion was considered 1.5 Pa, which corresponds to an erosion rate of 0.05 kg/m^{2} /s. Maximum erosion rate occurs at a distance of about 168 m from the shoreline (breaker line is about 172 m from shoreline). Certainly, it is expected that maximum erosion rate will occur near the breaker line, because, if the wave breaks, much of the wave energy will be dissipated and lead to significant erosion in the bed [12] .

Figure 15. Variation of erosion rate (less than 0.05 kg/m^{2} /s is surface erosion and more than 0.05 kg/m^{2} /s is mass erosion) versus distance to shoreline.

Figure 16 shows variation of the erosion rate versus bed shear stress. The low shear stress regime represents surface erosion, while the high shear stress regime corresponds to mass erosion. The top layer of loosely consolidated material overlies a layer of more consolidated material and can be easily disturbed and resuspended into the water column. For shear stress less than 0.5 Pa, erosion rate is low, while above 0.5 Pa, erosion rate increases rapidly with higher shear stress.

Figure 16. Variation of erosion rate versus bed shear stress.

The purpose of the current research was a better estimation of wave height and suspended cohesive sediment concentration profiles by a newly hybrid numerical model, which has been realized. Reviewing various applications of the present numerical model, different levels of accuracy in the prediction of wave height and concentration profile for cohesive sediment were observed. The models, such as Sanford and Maa’s erosion model (2001) [2] , Krone’s sedimentation model (1962) [3] , Tajima’s wave breaking model (2004) [4] , and the viscoelastic–plastic rheological model, were used to simulate interaction of the wave and muddy bed, and prediction of the concentration profile for cohesive sediment. The advection–diffusion equation was solved by the standard split method. Jazayeri Shooshtari and Soltanpour’s model and their reported experimental data, Sorourian’s model, Kemp’s field data and Sheremet’s field data were used for validation of the concentration profiles. Based on this research, the following specific conclusions can be drawn:
The generated effect by the presence of contaminants in the mud sediment is not included in the present model. A research on the combined behavior of fluid mud, transport and dissipation of contaminants is suggested as future study. In this study, cohesive sediments are used for concentration modeling, but it is recommended that a combination of fine and coarse aggregates is used for future studies.
Published on 06/10/16
Licence: Other
Are you one of the authors of this document?