Numerical modeling of the undertow structure and sandbar migration in the surfzone

A process-based 3D numerical model for surfzone hydrodynamics and beach evolution was established. Comparisons between the experimental data and model results proved that the model could effectively describe the hydrodynamics, sediment transport feature and sandbar migration process in the surfzone with satisfactory precision. A series of numerical simulations on the wave breaking and shoaling up to a barred beach were carried out based on the model system. Analyzed from the model results, the wave-induced current system in the surfzone consists of two major processes, which are the phase-averaged undertow caused by wave breaking and the net drift caused by both of the nonlinear wave motion and surface roller effect. When storm waves come to the barred beach, the strong offshore undertow along the beach suppresses the onshore net drift, making the initial sandbar migrate to the seaside. Under the condition of calm wave environment, both the undertow and net drift flow to the shoreline at the offshore side of the sandbar, and then push the initial sandbar to the shoreline. The consideration of surface roller has significant impact on the modeling results of the sandbar migration. As the roller transfer rate increases, the sandbar moves onshore especially under the storm wave condition.


Introduction
Sediment transport processes of the nearshore beach are directly responsible for the shoreline stability and coastal morphology evolution. Under the impact of global warming effect, recently, the occurrence of storms has greatly increased, bringing in more risk of potential beach erosion. In order to prevent the further loss of the valuable coastal resources, the development of an effective numerical modeling system for beach morphology evolution is necessary, and this topic has become increasingly popular these years.
The nearshore wave in the surfzone is an important process that directly affects the beach morphology. Once resuspended by waves, the wave-induced currents in the surfzone (e.g. the longshore current, undertow, and rip current) could transport the entrained sediments by patterns of both suspended load and bed load. Therefore, for obtaining reasonable modeling results of the beach evolution, the accurate prediction of the surfzone hydrodynamics is of great importance. The vertical structure of wave-induced undertow in the surfzone plays an important role for beach profile evolution since it directly controls the cross-shore sediment transport. Once the undertow develops, the re-suspended bed sediments by waves would transport either offshore or onshore, and then change the beach profile. Besides, the renewed beach profile also feeds back to the wave propagation and hydrodynamics in the surfzone. Generally, the sandbar migration is a result that caused by the mutually coupled processes among the surfzone hydrodynamics, sediment transport and morphology response.
Numerical modeling study of the surfzone hydro-sedimentological process is a classic and challenging topic in coastal researches. In previous works, the study methodologies mainly included the site observation and laboratory experiments, but the application of a general, process-based model was not many. The reason is that the traditional researches mainly focused on the 2DH models (e.g. Nam et al., 2011;Nam and Larson, 2010;Ding and Wang, 2008;Roelvink et al., 2009). To simulate the 2D wave-induced current structure, the concept of radiation stresses is widely applied as the driving force, but the concept of radiation stresses is under the depth-averaged system, which cannot present the vertical structure of the undertow. Although some researchers discussed the application of 3D models (e.g. Lesser et al., 2004;Roelvink, 2006), but in their models the driven force of the sediment transport is still tidal current, which could not model the vertical distribution of the wave-induced undertow structure. As a result, the present morphodynamic models cannot be directly applied for the modeling of beach profile evolution and sandbar migration.
Since the vertical distribution of the radiation stresses is a key factor which controls the development of the undertow in these years a number of formulae have been proposed to describe the depth-varying radiation stress, including the generalized Lagrangian mean (GLM) method (Ardhuin et al., 2008;Lin and Zhang, 2004), vertical mapping method (Mellor, 2005), vortices force method (McWilliams et al., 2004) and Eulerian mean method (Xia et al., 2004). Some researchers incorporated the above formulae into 3D hydrodynamic models to model the structure of the waveinduced current system (e.g. Wu and Zhang, 2009;Xie, 2011Xie, , 2012. Based on the ROMS platform, Warner et al. (2008) combined the 3D wave-induced current model and sediment transport model to simulate a typical tidal-inlet morphology evolution process, but their model did not focus on the cross-shore sediment transport, and the wave turbulent mixing effect was not considered in the model.
Generally, at present the application of a fully processbased 3D morphodynamics model system for the undertow and sandbar evolution study was still not reported. Because of this point, based on the platform of 3D wave-induced current model proposed by Xie (2011Xie ( , 2012, in this paper the sediment transport and morphology evolution modules were added to cast the model into a general, process-based surfzone morphodynamic model, and validations were carried out against the sandbar migration experimental cases (Arcilla et al., 1993). Finally, some discussions were given to further clarify the model performance.

Hydrodynamic model
The 3D numerical model of the surfzone hydrodynamics was developed on the platform of ECOMSED program (Xie, 2011). The governing equations are shown in Eqs.
(1)-(4). For modeling the wave-induced currents, the depthvarying wave radiation stresses, surface roller caused by wave breaking, and wave turbulent mixing effects are incorporated as the major driving forces. Cartesian coordinates are used in the horizontal and the terrain-following coordinate is used in the vertical.
where the vertical coordinate is the relative depth ranging from at the bottom to at the wave-averaged free surface; D=h+η is the total water depth; h is the water depth below still water level to the seabed; z is the vertical coordinate positive upwards with z=0 at the still water level; η is the wave-averaged free surface; t is time; x and y are the horizontal coordinates, respectively; U and V are the velocity components for x and y directions, respectively; ω is the velocity component under σ coordinate; g is gravity; p is pressure; M is the depth-varying wave radiation stress; R is the depth-varying surface roller momentum; K Mc and A Mc are the vertical and horizontal mixing coefficients combining waves and currents, respectively; ρ is the water density.
The formula proposed by Lin and Zhang (2004) is ap-plied for the vertical distribution of wave radiation stress, see Eq. (5).
where E=ρgH 2 /8 is the wave energy; n is the wave energy transfer rate; k is the wave number; T is the wave period; δ is the Kronecker symbol; i and j represent the x and y direction, respectively. The wave-current combined bottom shear stress is determined by Soulsby et al. (1993). The combined turbulent mixing coefficients are calculated by adding the current and the wave turbulent mixing coefficients (Xie, 2011(Xie, , 2012, see Eq. (6) and Eq. (7).
where A and K represent the horizontal and vertical turbulent mixing coefficients, respectively, and the subscripts M and W represent current and waves, respectively. The horizontal mixing coefficient A M for currents only is given by Smagorinsky formula, and the wave-induced horizontal turbulent mixing coefficient yields, where λ is the coefficient with a suggested range of 0.2-0.5. The current-induced vertical mixing coefficient K M is solved using a Mellor-Yamada 2.5 order turbulent closure model, and the wave-induced vertical mixing coefficient K W yields, where b is the coefficient with a suggested value of 0.001. Xie (2011) proposed an equation of the roller energy evolution, which yields, where α is the roller energy transfer rate; n is the wave direction vector; K R is the bed slope factor; ρ R and A R are the roller density and area, respectively; C r is the wave celerity; and C g is the wave group celerity.

R zn
The vertical profile of the roller momentum is expressed as Eq. (11). Because the depth integral of R z should be unity, it should be pre-normalized to following Eq. (12).
is solved, the roller energy can be calculated explicitly, and the corresponding terms in the governing equations caused by the roller can be determined, see Eqs. (13)-(15).
2.2 Wave model The combined refraction/diffraction wave model (REF/DIF) is used as the wave driver. It solves the parabolic mild-slope equation, and involves processes e.g. shoaling, refraction, energy dissipation, and irregular bottom bathymetry. The wave-current interaction can also be considered.

Sediment transport model
For the 3D non-cohesive sediment transport modeling system, the van Rijn (1984) theory was applied, which has been widely tested according to various literatures e.g. Lesser et al. (2004), Roelvink (2006) and Roelvink et al. (2009).
The governing equation of the suspended sediment transport model is shown in Eq. (16), where C is the suspended sediment concentration, D H and D V are the horizontal and vertical diffusion coefficients, ω s is the sediment settling velocity calculated by Eq. (17), D 50 is the medium sediment grain size, ν is water dynamic viscosity, ρ w is the seawater density, and ρ s is the sediment particle density.
. (17) The suspended sediments in the water could change the seawater density, and then bring in the density-driven baroclinic effect. In the model the seawater density should be updated in every time step following the expression in Eq. (18), where C vol is the sediment volume density. .
(18) The exchange between suspended sediments and seabed follows the concept of reference height, which is expressed as a=max (0.01D, 2.5D 50 ). The sediment concentration at the reference point is calculated using Eq. (19), where C a is the concentration at the reference point, T a and D* can be calculated following Eqs. (20) and (21), where τ cr is the critical shear stress for sediment entrainment, which can be calculated by Eqs. (22) and (23). The erosion and deposition fluxes E ro and D ep can be calculated following Eqs. (28) and (29), where kmx represents the center of the first vertical grid above the reference height.
u ′ * Below the reference point to the seabed, the sediment transport pattern is regarded as the bed load. Under the conditions of pure current motion and wave motion, the bed load transport rates can be expressed by Eqs. (30)-(33), respectively, where q bc and q bw are the bed load transport rates by current and waves, respectively, is the friction velocity, ub is the velocity at the first grid center from seabed, v r is the depth-averaged velocity corresponding to τ cr , U w is the maximum wave orbital velocity at the seabed. Based on the model results of the sediment transport rates, the bed elevation change could be expressed as Eq. (34), where z bed is the bed elevation.
2.4 Net drift caused by wave nonlinearity and roller effect Apart from the depth-varying radiation stresses, the strong nonlinearity of the shallow water wave also contributes to the wave-induced current system in the surfzone. As the incident waves arriving at the nearshore area, the action of the Stokes drift in the surfzone can bring in net movement of the suspended sediments over a wave period. Therefore, in the established model system the Stokes drift should be considered. Mellor (2003Mellor ( , 2005 derived the ana-lytical formulae of the vertical distribution of the Stokes velocity. Except the Stokes drift, the effect of the surface roller caused by the wave breaking process would bring in an additional residual flow, which also has potential impact on the suspended sediment transport. Based on the Mellor (2005) expressions, Warner et al. (2008) and Haas and Warner (2009) added the roller effect, see Eqs. (35) and (36). Note that Eqs. (35) and (36) contain both the Stokes drift and the roller-induced residual flow, thus in this paper the combined velocities are called as "net drift".
In the original expressions of Eqs. (1)-(4), the driven forces of the wave-induced current motion do not include net drift, so Eq. (37) should be applied to consider the impact of net drift, where the superscripts M, O, and N represent the modified, original and net drift velocities, respectively.
(37) When solving the governing equations, the original velocities are used in the hydrodynamics model, while the modified velocities are used in the sediment transport model as well as the bed shear stresses. This method has been widely applied in Warner et al. (2008), Haas and Warner (2009) and Roelvink et al. (2006).
To better clarify the model procedure, Fig. 1 illustrates the flow chart of the model system. Once the simulation begins, the wave module is firstly operated, and it provides the depth-dependent wave radiation stresses, surface roller momentum, net drift velocity and wave mixing coefficient to the hydrodynamic model. After solving the hydrodynamics, the velocity and wave-current combined bottom shear stresses are determined, and then the sediment transport and morphology modules begin to run. For the next time step, the bathymetry is updated based upon the bed level change, and the density term is also renewed according to the new suspended sediment concentration.

Validation of the numerical model
In order to test the behavior of the established model system, the LIP11D experiment Cases 1b and 1c were applied (Arcilla et al., 1993). The experiments were carried out in the large wave flume in Delft Hydraulic Institute. The length, width and height of the wave flume are 233 m, 5 m and 7 m, respectively.
As to the experimental conditions, in Case 1b the strong waves has been acted on the initial barred beach for 18 h with the incident wave height of 1.4 m. Based on Case 1b, the beach erosion and the offshore migration of the sandbar were investigated, which represented the beach profile evolution under a typical storm process. In Case 1c, the weak waves of 13 h were generated with an incident wave height of 0.6 m. Based on Case 1c, the onshore sandbar migration was investigated, which represented the beach profile recovery under calm weather. In both cases (Case 1b and Case 1c), the mean sediment grain size D 50 is 0.22 mm. Table 1 lists the main parameters in the experiment. Fig. 2 presents the experimental bathymetry of the initial sandbar as well as the location of the observation stations. Table 2 lists the parameters applied in the numerical model.
The experimental environment of Case 1b represents the storm wave condition, and the incident wave height at the beach foot reaches 1.6 m, thus for this case the waves continuously break along the flume, and the location of the wave breaking point is near x=40 m. Owing to the strong wave breaking, the maximum wave setup is approximately 0.06 m at the shoreline. The experimental environment of Case 1c represents the calm wave condition, and the incident wave height at the beach foot is only 0.6 m, thus the wave breaking effect is much weaker than that in Case 1b, the location of breaking point in case 1c occurs until x=130 m (near the sandbar top), and the maximum wave setup is only 0.03 m at the shoreline. For both cases, the wave height slightly increases in the trough between the two initial sandbars (140 m<x<160 m). Fig. 3 and Fig. 4 present the comparisons between the measured and simulated wave height and wave setup values at t=1 h, respectively. According to the comparisons, both the modeled wave height and wave setup are in satisfactory agreement with the observed data, which indicates that the selected wave solver has a good performance for the cases applied in this paper. Fig. 5 and Fig. 6 give the comparisons between the simulated and measured phase-averaged current speeds for Case 1b and Case 1c at t=1 h, respectively. Based on the comparisons, the simulated current speed values are in good agreement with the measured data, showing that the model can effectively describe the vertical structure of the undertow at different locations.
For Case 1b, the surfzone width is broad due to the strong breaking process, and the surface roller effect is more    significant. Therefore, in Case 1b the maximum returning offshore undertow speed could exceed 0.3 m/s especially over the sandbar top, and near the seabed the undertow direction points offshore along the flume. For Case 1c, the incident waves are weak, thus the undertow speed is smaller than 0.2 m/s. Different from Case 1b, at the offshore side of the sandbar, the undertow direction near the seabed points to the shoreline. The detailed discussion of such differences between the two cases is given in Section 4.1. Fig. 7 presents the comparison between the measured suspended sediment concentration and the simulated results for Case 1b at t=1 h. The results show that the agreement between the two datasets is good in the offshore area where the beach profile is smooth (60 m<x<130 m), and the maximum errors occur at the offshore side of the initial sandbar (130 m<x<140 m). Around the sandbar, the model slightly underestimates the sediment concentration values. Actually, in the experiment many small sand ripples might develop in the surfzone especially around the sandbar, which could make the beach profile not smooth. The existence of sand ripples at the seabed would generate small-scaled turbulence eddies near the ripples, and may entrain the bed sedi-  ments around the ripple, but the phase-averaged model cannot account for this sub-grid process. Besides, the bed roughness height was set to be constant along the flume (2.5D 50 ) in the numerical model, so the roughness variation caused by such ripples was neglected. Fig. 8 gives the comparisons between the measured and simulated beach profile and sandbar location for Case 1b (at t=18 h) and Case 1c (at t=13 h), respectively. The model is shown to well describe the offshore and onshore migration process of the sandbar for both cases.
Generally, combining the above comparisons, it proves that the model can effectively present the undertow structure and the offshore/onshore sandbar migration processes in the surfzone with a satisfactory precision.

Discussions on the hydro-sedimentological process
In the above section, the behavior of the established model has been validated using the LIP11D experimental datasets. In order to shed more lights on the model performance, in this section some discussions are given based on the model results, and a sensitivity analysis of the key calibration parameter of the model is also made.

Distribution of the undertow and sediment transport
In the LIP11D experiments, the conditions of case1b and case 1c respectively represent the offshore and onshore migration processes of the initial sandbar. Fig. 9a and Fig.  9b compare the distribution differences of the modeled undertow fields between the two cases (t=1 h for example).   XIE Ming-xiao et al. China Ocean Eng., 2017, Vol. 31, No. 5, P. 549-558 555 Some discussions are made as follows.

∂η/∂x
(1) Case 1b presents the offshore sandbar migration under storm waves, and the width of breaking zone nearly covers the whole beach profile owing to the strong incident waves. According to the experiment data and simulation results, the wave setup continuously increases (see Fig. 4a), and that means the direction of the mean water level pressure gradient is approximately uni-directional along the beach. As a result, along the whole beach profile, the water in the surface layers moves onshore to balance the water level pressure gradient, and consequently forms a vertical returning undertow flowing offshore near the bottom, which means that only a single phase-averaged vertical circulation develops and covers all the surfzone. Near the top of the initial sandbar, the offshore undertow speed reaches the maximum value due to the strongest breaking process. Therefore, under the continuous action of the strong offshore undertow, the sediments suspended by breaking waves at the sandbar top can be carried and transported to the seaside, and consequently lead to the offshore sandbar migration.

∂η/∂x
(2) Case 1c presents the onshore sandbar migration under calm wave environment. For this case, the surfzone width is much narrower than that of case 1b. Based on the modeling result shown in Fig. 4b, the mean water level presents the feature of setdown at the seaside of the sandbar, and then turns to setup owing to the wave breaking at the sandbar top where the water depth quickly decreases. Therefore, the value of mean water level pressure gradient is opposite at the offshore and onshore sides of the initial sandbar, which then makes the undertow structure totally different at each side of the sandbar. At the seaside of the sandbar, the phase-averaged undertow direction near the seabed points to the shoreline, and the undertow direction at the surface points offshore, but at the onshore side of the sandbar, the undertow direction near the seabed turns to point offshore, and the undertow direction at the surface points to the shoreline. This feature indicates two opposite phase-averaged vertical circulations near the sandbar, which is different from the uni-directional undertow structure of Case 1b.
(3) Apart from the returning undertow driven by the mean water level pressure gradient, the wave period-averaged net drift also contributes to the cross-shore sediment transport. The direction of the net drift is onshore, and the drift velocity gradually decreases from the surface to the seabed. Therefore, the phase-averaged flow direction along the beach is actually determined by the balance between the undertow and net drift. For Case 1b, the strong offshore undertow near the seabed suppresses the onshore net drift, so the net flow direction points offshore along most part of the beach. But for Case 1c, the phase-averaged flow directions outside the surfzone for both undertow and net drift are the same, then the onshore velocity near the sandbar is en-hanced, which makes the sandbar move onshore.
(4) Fig. 10 compares the differences between the total sediment transport rate along the beach for Case 1b and Case 1c at t=1 h. According to the above analysis of the hydrodynamic features, it indicates that under the continuous wave breaking process of Case 1b, the sediment transport direction is mainly offshore, while for Case 1c the sediment transport direction is mainly onshore. For both cases, the maximum transport rate value appears near the sandbar top. Actually, once the bed level elevation changes, the sediment transport rate along the beach profile is also updated.
Generally, the sediment transport and sandbar migration processes in the surfzone are very complex, and they are directly determined by the vertical structures of both undertow and net drift. According to the modeling results, for the wave-induced current system in the surfzone, there are two different mechanisms which are: (1) the returning undertow near the seabed caused by the mean water level pressure gradient, and (2) net drift caused by both the wave nonlinearity and roller residual flow. Under the condition of storm wave case, the wave breaking is fully developed, which induces the strong offshore undertow near the seabed, and the onshore net drift is suppressed, which then results in the offshore sediment transport. While under the condition of calm wave case, the phase-averaged undertow direction points to the shoreline at the seaside of the sandbar, and then pushes the bed sediments to the shoreline. Actually, in the realistic nearshore environment, these two processes (undertow and net drift) coexist and mutually interact, and make the sandbar migration characteristics extremely complicated. Therefore, in the modeling of the sediment transport process in the surfzone, the depth-dependent radiation stresses, surface roller effect, and net drift should be taken into account simultaneously.
Based upon the modeling results and analysis, the mechanism of the wave-induced currents and cross-shore sediment transport characteristics are discussed. Additionally, some researchers also proposed the streaming process in the wave boundary layer should be another potential hydrodynamic factor leading to the onshore sediment transport at the seabed (e.g. Longuet-Higgins, 2005;Tajima and Madsen, 2006). The model established in this paper is under the phase-averaged system, which does not account for the sub- wave process, so the streaming effect could not be easily considered. At present, the accurate modeling of the streaming process and its impact to the bed sediment transport is still a challenge for coastal modelers. However, from the comparison with the observed datasets, the model could provide the basic hydrodynamic background in the surfzone, and the model results are in good agreement with the measured data.

Sensitivity analysis of the surface roller effect
After the establishment of a numerical model system, the judgmental parameters in the model, especially the key calibration parameters should be analyzed. For the non-cohesive sediment transport model, once the sediment pattern is given, the related physical parameters e.g. the mean grain size, settling velocity, critical shear stress for erosion are then determined and should not be simply calibrated. Therefore, for the specified cases of sandbar migration applied in this paper, the accurate hydrodynamics background is the most important.
For the key calibration parameters in the hydrodynamics model, Xie (2011) has given a detailed discussion, and suggested their values. Since the cases used in this paper are under the 2DH frame, so the key judgmental parameters which control the vertical distribution of the undertow are (1) coefficient b used in Eq. (12), and (2) the surface roller transfer rate α and roller density ρ R used in Eq. (13). Through a series of sensitivity analysis, b is suggested as 0.001 after validation using a large series of experimental datasets, and the roller transfer rate is regarded far more sensitive than the roller density (Xie, 2011).
Actually, the surface roller effect accompanies with the wave breaking process, so it is a special and important process in the surfzone, and the momentum transfer of the roller energy dissipation directly contributes to the development of the undertow as well as the roller-induced residual flow. Therefore, the roller transfer rate α should be the most sensitive and judgmental parameter in the model calibration. From this aspect, in this section the sensitivity analysis of roller transfer rate is mainly focused on.
In the sensitivity analysis, two characteristic values of the roller transfer rate α were added, and they are α=0.0 which represents the effect of no roller, and α=1.0 which represents that the roller is fully developed. Fig. 11a and Fig. 11b illustrate the simulation results with different α values for Case 1b and Case 1c, respectively. Some discussions are given as follows.
(1) For Case 1b, the storm waves induce a significant breaking process nearly covering the entire barred beach, thus the surface roller is developed simultaneously with the breaking waves. From the comparison shown in Fig. 11a, the value of the roller transfer rate α has a significant impact on the modeled beach profile. As α increases, the top of the sandbar further moves to the shoreline. When the roller is fully developed, the location of the maximum undertow speed also moves to the shoreline, and then further changes the sandbar evolution. Based on the sensitivity analysis, the roller transfer rate α=0.45 is the proper value, and this value is in good agreement with that proposed by Xie (2011), which suggested that α should be in the range of 0.2-0.5.
(2) For Case 1c, generally the impact form of the roller transfer rate is the same as that of Case 1b. However, because the intensity of wave breaking for Case 1c is much weaker than that of Case 1b, thus the degree of impact is also smaller than that in Case 1b. Actually, the mechanism of the surface roller usually accompanies with the wave breaking, therefore, as the breaking becomes more intense, the roller effect would then increase.
Based on the above analysis, in the numerical modeling of the sandbar migration, the surface roller effect should be considered especially under the storm case, and a reasonable value of the roller transfer rate should be carefully calibrated.
To sum up, in this paper we have established an integrated modeling system for the surfzone hydro-sedimentological processes, and the validation results indicate that the model has satisfactory behavior. However, in the present work the Lesser et al. (2004) theory is applied as the sediment transport module, and in future works we will further focus on the comparison between other theories.

Conclusions
(1) A process-based 3D surfzone hydrodynamics and beach morphology evolution model was established, in which the processes of wave-induced undertow, suspended and bed sediment transport as well as the morphology re- sponse were incorporated. A series of wave flume experimental cases were applied to validate the model. The results proved that the model could describe the undertow structure and the sandbar migration characteristics in the surfzone with a satisfactory precision.
(2) The vertical wave-induced current system in the surfzone consists of two main processes, which are the phase-averaged returning undertow caused by wave breaking and net drift caused by both the wave nonlinearity and roller-induced residual flow. When storm waves come to the barred beach, the strong offshore undertow along the beach suppresses the onshore net drift, and makes the initial sandbar migrate to the seaside. Under the action of calm waves, both the undertow and net drift point to the shoreline at the offshore side of the sandbar, and then push the initial sandbar to the shoreline.
(3) The surface roller accompanies with the wave breaking, and its evolution is a special and important process in the surfzone. The value of roller transfer α has a significant impact on the beach profile evolution in the numerical model. As α increases, the modeled sandbar moves to the shoreline, and this effect becomes even obvious when the wave breaking is more intense. Therefore, in the modeling of the sandbar migration, the surface roller effect should be taken into account. This is especially true for the storm case. Therefore, the reasonable value of the roller transfer rate should be carefully calibrated in the modeling process.