Investigation and Discussion on the Beach Morphodynamic Response Under Storm Events Based on A Three-Dimensional Numerical Model

A well-established 3D phase-averaged beach morphodynamic model was applied to investigate the morphodynamics of a typical artificial beach, and a series of discussions were made on the surfzone hydro-sedimentological processes under calm and storm events. Model results revealed that the nearshore wave-induced current presents a significant 3D structure under stormy waves, where the undertow and longshore currents exist simultaneously, forming a spirallike circulation system in the surfzone. Continuous longshore sediment transport would shorten the sediment supply in the cross-shore direction, subsequently suppress the formation of sandbars, showing that a typical recovery profile under calm waves does not necessarily develop, but with a competing process of onshore drift, undertow and longshore currents. Sediment transport rate during storms reaches several hundreds of times as those under calm waves, and two storm events contribute approximately 60% to the beach erosion. Sediment transport pattern under calm waves is mainly bed load, but as the fine sands underneath begin to expose, the contribution of suspended load becomes significant.


Introduction
Affected by global warming and climate change, the occurrence of storm surges increases in recent years, bringing in additional risks of the nearshore beach erosion. To prevent further loss of the coastal resources, the development of an effective numerical modeling system for beach morphology evolution is mandatory. Situated in the surfzone, the strong wave nonlinearity and breaking make the hydrodynamic process extremely complex, and the wave velocity asymmetry leads to the non-equilibrium sediment transport. Besides, the wave-induced current system in the surfzone presents various features e.g. the longshore current, rip current, undertow, surface roller, Stokes drift and boundary layer streaming, and these mutually coupled processes form a 3D phase-averaged current system (Svendsen and Lorenz, 1989). As a result, for a typical morphodynamic model, the accurate description of the 3D wave-induced current structure directly determines the model precision and its validity.
However, over the last decades, most of the model systems focused on the 2DH models (Zhang et al., 2019;Nam et al., 2011;Nam and Larson, 2010;Roelvink et al., 2009;Ding and Wang, 2008). The concept of radiation stresses initially proposed by Longuet-Higgins (2005) was incorporated as main driving force of the wave-induced currents, but the radiation stresses are under the depth-averaged frame, which cannot present the vertical structure of the undertow. In these years, a number of formulations have been proposed to describe the depth-varying radiation stresses e.g. (Mellor, 2005), vortices force method (McWilliams et al., 2004) and Eulerian mean method (Xia et al., 2004). Some researchers incorporated the above formulations into 3D hydrodynamic models to simulate the structure of the wave-induced current system (Wu and Zhang, 2009;Xie, 2011Xie, , 2012. In the past decade, the combination with the morphodynamic model becomes popular, e.g. Nam et al. (2010) proposed a 2DH model to simulate the tombolo evolution, Roelvink et al. (2006) developed the XBeach package. Based on different platforms of DELFT3D-FLOW and ROMS, Lesser et al. (2004) and Warner et al. (2008) have respectively developed the 3D morphodynamic models, but their models did not include a sophisticated 3D wave-induced current module. Xie (2011) has originally developed a fully 3D wave-induced current model, in which several key processes including the depth-varying radiation stress, net drift, surface roller evolution and turbulent mixing effect were incorporated. The above model has been experimentally verified to be well suited for presenting the characteristics of various nearshore processes e.g. the wave setup/down, undertow structure, and longshore currents. Xie (2012) further tested the model behavior under a 3D case including rip currents and meandering longshore currents. By incorporating sediment transport and morphology evolution modules, Xie et al. (2017) casted the model system into a generalized surfzone morphodynamic model. Validated by the datasets from an experiment carried out in the large wave flume (Arcilla et al., 1993), the model was also shown to be effective to describe the sediment concentration profile in the surfzone, as well as the onshore/offshore sandbar migration beneath both stormy and calm waves.
Nevertheless, considering that the works presented by Xie et al. (2017) only covered the 2DV case of beach profile evolution in a flume, in this paper a real-life beach engineering application was applied to further evaluate the model behavior under the full-scale site case. Firstly, a series of simulations have been run to repeat the hydro-sedimentological process as well as the beach evolution, which can be regarded as an extended validation of the established model system. After that, a series of discussions are expanded according to the model results, including the 3-D structure of surfzone hydrodynamics, as well as the sediment transport features.

Hydrodynamic model
The 3D numerical model of the surfzone hydrodynamics was initially developed by Xie (2011). The governing equations are shown in Eqs. (1)−(4), amongst which the depth-varying radiation stresses, surface roller caused by wave breaking, and wave turbulent mixing effects are incorporated as major driving forces. (2) where the vertical coordinate is relative depth; is total water depth; h is water depth below still water level to seabed; z is vertical coordinate positive upwards with z=0 at still water level; is wave-averaged free surface; t is time; x and y are horizontal coordinates, respectively; U and V are velocity components for x and y directions, respectively; is velocity component under coordinate; g is gravity; p is pressure; M is depth-varying wave radiation stress; R is depth-varying surface roller term; K Mc and A Mc are vertical and horizontal mixing coefficients combining waves and currents, respectively; is water density. The formulation proposed by Lin and Zhang (2004) is applied for the vertical distribution of wave radiation stress, see Eq. (5), where is wave energy; n is the vector of wave energy transfer rate; k is wave number, with k x and k y are components in x and y directions, respectively; T is wave period; is Kronecker symbol, where the subscripts i and j represent x, y direction, respectively.
where A and K represent 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 formulation, and the wave-induced horizontal turbulent mixing coefficient yields, λ where is a coefficient with a suggested range of 0.2−0.5. The current-induced vertical mixing coefficient K M is solved by using a Mellor-Yamada 2.5 order turbulent closure model, and the wave-induced vertical mixing coefficient K W yields, where b is a coefficient with a suggested value of 0.001. As mentioned above, besides the radiation stresses, an aerated "white capping" over the top of breakers could bring in additional momentum transfer to the surfzone hydrodynamic processes, which is recognized as surface roller. Xie (2011) derived an energy balance equation of such roller, which yields, where C r is wave celerity, C g is wave group celerity, is the roller energy transfer rate, K R is slope factor, is roller density, A R is roller area. 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 R zn following Eq. (12). Once A R 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-induced net drifts As waves approaching the nearshore area, it is well-acknowledged that the Stokes drift can bring in net mass transport over a wave period. Therefore, in the established model the Stokes drift should be considered. Mellor (2003Mellor ( , 2005 proposed an analytical formula of the vertical distribution of the Stokes velocity. Additionally, the effect of the surface roller caused by the wave breaking process would bring in an additional residual flow (roller drift), which also has potential impact on the mass transport. Warner et al. (2008) and Haas and Warner (2009) added such effect, see Eqs. (16)−(17). Note that Eqs. (16)−(17) contain both the Stokes drift and the roller-induced residual flow, thus in this paper the combined velocities are called as "net drifts".
In the original expressions of Eqs.
(1)−(4), the driving forces of the wave-induced current motion do not include the net drifts, as a result Eq. (18) should be applied to take account their contribution, where the superscripts M, O, and N represent modified, original and net drift velocities, respectively.
(18) 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, known as the GLM treatment, has been applied in pre-existent studies e.g. Warner et al. (2008), Haas and Warner (2009) and the XBeach model (Roelvink et al., 2006).
Besides the undertow and net drift in the upper waters, another potential wave-induced process is the steady streaming within the bottom boundary layer which may contribute to the onshore sediment transport. Note that the non-equilibrium sediment transport within the boundary layer is not easy to incorporate in a typical phase-averaged model, the treatment of steady streaming follows that proposed by Walstra et al. (2000), in which the streaming is modeled as additional shear stress.

Wave model
In this paper, we dedicate to developing a phase-averaged hydro-sedimentological model. Under this case, the parameterized wave model seems more suitable for providing the wave information. Therefore, we applied the wellknown SWAN model as the wave driver, which involves several key processes e.g. shoaling, refraction, diffraction and energy dissipation.

Sediment transport model
ν The governing equation of the suspended sediment transport model is shown in Eq. (19), where C is suspended sediment concentration, D H and D V are horizontal and vertical diffusion coefficients, respectively. There is a large quantity of formulations in the literature for calculating sediment settling velocity, the expression of which used in this paper is formulated by Wang (2001), see Eq. (20), where D 50 is the medium sediment grain size, and is water dy-namic viscosity.
The suspended sediments in the water could change the seawater density, and then bring in the potential densitydriven baroclinic effect. In the model the water density updates in every time step following the expression shown in Eq. (21), where C vol is the sediment volume density, and are seawater and sediment particle densities, respectively. .
τ cr τ cw 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 by Eq. (22), where C a is the concentration at the reference point, T a can be calculated by Eq. (23).
is the critical shear stress for sediment entrainment, and is the wave-current combined bottom shear stress as formulated by Soulsby et al. (1993) and Xie et al. (2021). The erosion and deposition fluxes E ro and D ep can be calculated by Eqs. (24)−(25), where kmx represents the center of first vertical grid above the reference height. .
(25) Below the reference point to the seabed, the sediment transport pattern is regarded as bed load. Under the conditions of pure current motion and wave motion, the bed load transport rates can be expressed by Eqs. (26)−(29), respect-V ively, where q bc and q bw are bed load transport rates by current only and wave only, respectively, u * is the friction velocity, u b is the velocity at the first grid center from seabed, is the depth-averaged velocity, 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. (30), where z bed is the bed elevation, Δz is the vertical distance from the reference level a to the grid point nearest to the bed.
Having developed from the open source ECOMSED platform, the surface boundary conditions, including velocity, turbulent and sediment concentration in the model follow the standard zero-gradient treatment, and at the bottom a no-slip boundary condition is imposed. Readers of interest could refer to the ECOMSED user manual for more details. As for the sediment transport model, the upper boundary condition follows the zero-gradient treatment, and at the reference level near the bottom, let C = C a .

Morphology evolution characteristics of the beach
A real-life engineering case, the Weifang artificial beach was selected to test and discuss the model performance. The artificial beach situates at the bottom of the Laizhou Bay, which is a semi-closed bay in the Bohai Sea, China. The geophysical location of the studied beach is presented in Fig. 1, and Fig. 2 shows its layout. The artificial beach was constructed in October 2011, with 3.3×10 5 m 3 of sea sands carpeted to create the beach. Two breakwaters of 250 m were constructed at each side of the beach, and the beach foot line presents a "wave-like" shape, which is broader at the western side. Beach slope ranges from 1:50−1:25 with top elevation of 4.0 m above MSL, and the bed elevation at the sea side of the beach is from −1 m to 2 m. The beach consists of two different sediment types vertically (coarse sand at the surface and fine sand underneath), at the surface layer the mean sediment grain size D 50 is 0.67 mm with the thickness of 0.5 m, and beneath the surface layer the D 50 is 0.13 mm. According to the field survey carried out by Hou et al. (2013) around the study site area, the offshore bed sediments have a D 50 of 0.12 mm.
Once finished, the artificial beach had faced severe erosion after the coming winter and spring seasons (from October 2011 to May 2012, nearly half a year), and the initially carpeted beach sands experienced a loss of approximated 67%. Hou et al. (2013) analyzed the data of the adjacent Weifang Harbor, and proposed that the astronomic tide presented a regular semi-diurnal feature with tidal range only 1.6 m (mean HWL 0.73 m and mean LWL −0.87 m). Using the FVCOM modeling package, Liu et al. (2013) simulated the tidal current motion in Laizhou Bay, and concluded that the mean tidal current speed near the artificial beach was only 0.11 m/s during spring tidal cycle, and even decreased to 0.06 m/s during neap tidal cycle. Commonly, for a beach case, the tidal current speed is far less than that of the wave-induced currents in the surfzone, which cannot directly re-suspend the beach sands. Therefore, the key dynamic factors leading to the beach morphology evolution should be the combined actions of the nearshore waves and wave-induced currents.
Tianjin Research Institute for Water Transport Engineering (TIWTE) has installed a temporary wave observation station at the water depth of 9.5 m (see Fig. 1b), and wave data from December 2008 to January 2010 were recorded. Fig. 3 illustrates the annual distribution of the wave frequency, from the figure we can see that the dominant wave directions are from NE and NNE with frequencies of 25.3% and 17.4%, respectively, and strong waves mainly come from the NE direction. Fig. 4 illustrates the measured beach elevation variation from October 2011 to May 2012. Once constructed, the entire beach has been significantly scoured especially at the middle and east parts with maximum elevation change exceeding 3.0 m. The erosion was less intense at the western end, and the beach slightly deposited at the very end of the western breakwater with thickness of about 1.0 m. Based on the evaluation by Zhou (2014), the volume of sediment loss of the beach is 22.1×10 4 m 3 , which takes a percentage of 67% compared with the initially carpeted sands. The tidal current around the beach is extremely weak, thus the key factors affecting the beach erosion should be the nearshore waves and wave-induced currents. When the incident waves become intense, the beach sediments can be re-suspended, and then be transported in both alongshore and cross-shore directions. Fig. 3 shows that the dominant offshore wave direction is from NE, which can directly approach the beach without any obstruction. The crest of the NE waves has a large angle with the shoreline, which could induce the continuous longshore sediment transport even under calm weather. As the beach is situated in northern China, in winter and spring the cold air coming from higher latitude region brings in storm surges, and subsequently has significant impact on the Weifang coast. Most important of all, the artificial beach is located at the bottom of the semi-closed Laizhou Bay, thus the water level setup caused by the storm surge could be enlarged even the winds are not particularly strong, see Table 1.

Modeling and validation of the beach morphology variation
4.1 Model setup Considering the existence of three different sediment patterns of the surface beach (0.67 mm), bottom beach (0.13 mm) and offshore bed sediment (0.12 mm), in the model we separate these sediments corresponding to their locations (horizontally and vertically). At the beginning, the sediments grain size is set to 0.67 mm at the beach and 0.12 mm elsewhere, and as the simulation proceeds, the sediment grain size changes to 0.13 mm when the beach erosion reaches a threshold of 0.5 m. It should be stressed that in the present model the mixing effect between these different sediments cannot be taken into account due to its complex nature.
According to the analysis of the beach erosion in Section 3, the offshore waves coming from NNE−NE play a dominant role in the beach evolution, and the tidal currents could be neglected for this specified case. The key parameters applied in the numerical model are presented in Table 2. Besides, the near-bed undertow should be particularly considered when modeling the sediment transport in the crossshore direction (Xie et al., 2017), thus the vertical layers near the seabed is refined.
The tidal water level applied in the simulation only varies with time throughout the model domain. There are two reasons why we follow this method: (1) Compared with the strong longshore currents caused by wave breaking, the tidal current is very weak, the mean speed of which is only 0.11 m/s during spring tidal cycle, and even decreases to 0.06 m/s during neap tidal cycle (Liu et al., 2013). Hence, the dominant flow pattern leading to the beach morphology evolution should be the wave-induced currents.
(2) The model domain only covers a small spatial range of 2400 m× 1200 m, and thus the tidal level has trivial variation throughout the model domain.
For the specified case, besides the surface and bottom conditions described in Section 2.4, in the consideration of lateral boundary conditions, we use the zero-gradient treatment for all variables, including the velocity, sediment concentration and sediment transport rate. That means no additional sediment supply is considered outside the model domain, which is reasonable because the studied beach is arti-ficial.

Representative condition
The environment of coastal waves has stochastic nature, and accordingly we cannot obtain the realistic wave information in every computation time step. Therefore, the concept of representative condition is mandatory for modeling the long-term morphodynamic processes. Steijn (1989) has proposed two methods, which are the MRF (multiple representative waves) and SRF (singular representative waves). According to the annual wave rose shown in Fig. 3, the offshore wave directions that contribute most to the beach erosion highly concentrate in the range of NE−NNE. Additionally, after refracting in the shallow water, the offshore waves incident from these two directions (NE and NNE) could get even closer at the beach foot. Therefore, the SRW method is favored in this case. Based on the wave direction spectrum, we weight the occurrences of these two major directions and have the representative wave direction of 39.6°. The selection of representative wave height not only relates to the wave spectrum, but also depends on which physical process is concerned. Steijn (1989) and Cayocca (2001) proposed that the longshore sediment transport is proportional to H 2.5 , and the cross-shore sediment transport is proportional to H 6 . In this paper, the erosion of the artificial beach is mainly caused by the continuous longshore current, thus we use H 2.5 as the principle to calculate the representative wave height, see Eq. (31), where f means the frequency, i means different wave height grade.
The representative wave period could be calculated following the wave height-period joint distribution function analyzed from the field measurement data, which is set to 4.7 s. Validated by the above wave parameters, a regional SWAN wave model was set up to provide the wave field near the artificial beach area, and the wave spectrum is selected as the JONSWAP pattern.
Storm surges have significant impacts on beach erosion, which must be considered thoroughly in the model. According to historical weather data, the Laizhou Bay has encountered several storm surges in the modeling duration. Sheltered by the adjacent long breakwater of Weifang harbor, the storm waves from NNE−ENE could directly affect the beach, which should gain special attention in the modeling process. We have found that the two cold wave weather processes are directly responsible for the beach erosion, which are (1) December 14 2011 to December 15 2011, and  Liu et al. (2013) modeled the two storm surge processes applying FVCOM and SWAN models, and their results could provide the water level and offshore wave parameters to the present beach morphodynamic model. Fig. 5 illustrates the water level and wave height processes at the seaward side of the beach foot (−5 m contour, see Fig. 2). According to the model results, during the storm surge the highest water level could reach 3.31 m, which agrees well with that analyzed by Zhou (2014). Table 3 gives the arrangement of representative conditions, in which periods C1−C3 mean calm waves and S1−S2 mean two storm surges, respectively.

Model results and discussions
Beach morphodynamics in the surfzone has an extremely complex nature, and especially several processes e.g. the bed level variation, wave field, wave-induced currents and sediment transport are strongly coupled. Therefore, once the model error occurs when describing a specified dynamic factor, the errors tend to be quickly accelerated, which makes the model no longer robust. Therefore, to test the model behavior, if we only check the results of the modeled beach elevation variation, it is not convenient for us to distinguish the "initial error".
Considering the above reason, we carried out modeling of two different cases, and they are the fixed bed case (Case A) and the movable bed case (Case B). In Case A, the sediment transport and bed evolution modules are switched-off, that is, only the hydrodynamics are modeled without the bed level feedback. The purpose of Case A is to evaluate the validity of the modeled wave field and wave-induced current system. Once the modeled hydrodynamics are recognized to be acceptable for driving the sediment transport module, in Case B we then model all hydro-sedimentological processes and finally obtain the beach morphology evolution. By combining the above two steps, we can better test the model performance.

Waves propagation
By running Case A, Fig. 6 presents the modeled wave fields for several water level and incident wave conditions. Under calm weather, when the water level is low, the waves break only at the beach foot, and most of the beach area is dry, which cannot lead to erosion. As the water level rises, the breaking line moves shoreward, leading to a larger erodible beach width. When facing the storm surges, the water level significantly increases under the impact of the setup, and nearly the entire beach could be submerged beneath the breakers (Fig. 6d), which can result in strong beach erosion.

Wave-induced currents
By running Case A, Fig. 7 presents 2DH wave-induced current fields under different water level and wave conditions. Strong currents highly concentrate in the nearshore region, and quickly diminishes offshore. With higher water level and stronger offshore incident waves, the surfzone gradually moves onto the beach. Under calm waves, at MSL the maximum longshore current speed ranges from 0.3 to 0.4 m/s, and at MHL the current speed ranges from 0.4 to 0.7 m/s. The longshore current has greater intensity from the middle to western part of the beach, which is because the wave incident angle is more oblique to the shoreline, and the eastern breakwater could block the offshore waves, which induces a reversed circulation at the eastern end of the beach.
The horizontal distribution of the longshore currents well explains the present beach erosion. Based on the model results, the depth-averaged wave-induced current speed reaches 0.7 m/s under calm weather, and even exceeds 1.5 m/s during the storm surge. Once entrained by the nearshore waves, the beach sands could be continuously transported westward by the strong longshore currents, which induces the sand loss.

Bed elevation variation
Proving that the model could provide reasonable hydrodynamic background, subsequently we carried out movable bed case (Case B), and Fig. 8 illustrates the beach elevation variations after two storm surges. The entire beach encounters significant erosion during storm surges. After storm surge #1, the maximum scour thickness of the middle and eastern parts of the beach (sections P5−P6) exceeds 2.5 m, meaning that the fine sediments of the lower beach fully expose to strong waves. Trapped by the western breakwater, the beach slightly deposits with a thickness of 0.5 m around section P1. In the offshore, the bed elevation raises about 1.0 m, indicating a seaward sediment transport by the undertow. The horizontal distributions of the erosion/deposition after both storm #1 and #2 are similar, implying that the beach evolution during the calm wave conditions does not play a crucial role in the present beach erosion. Fig. 9 presents the comparison between the modeled and measured beach profiles at different representative locations (P1−P8), and two model skill parameters e.g. the mean  XIE Ming-xiao et al. China Ocean Eng., 2021, Vol. 35, No. 1, P. 12-25 square deviation (RMS) and correlation coefficient (COR) are given to evaluate the model error. Generally, it shows that the modeled beach profiles agree well with that observed with COR=69%−84%, and RMS=0.32−0.89 m. In general, the modeled beach elevation has larger deviation at the eastern part of the beach. Besides, the model also underestimates the beach erosion at the western side, while overestimates the beach erosion near the eastern breakwater. The potential reasons that may induce the errors are discussed as follows.
(1) The beach morphology responses quickly to the incident waves as well as the nearshore currents. When lack-ing in the time-series measured wave data, the wave condition applied in the model has to be simplified. That is, only two strong storm surges from NE−ENE as well as the representative calm wave condition were taken into account. However, in addition to these two selected storm surges, during the modeled time span the Laizhou Bay had encountered several other cold wave weather processes with smaller intensity, as well as other wind directions e.g. from NW−N. Therefore, the inputted wave conditions in the model did not cover the whore information. Additionally, despite the offshore waves from N−NW which could be largely sheltered by the breakwater of the Weifang Harbor,  the locally generated wind waves still have impacts at the eastern part, which may generate the eastward longshore sediment transport, and that explains why the model overestimates the beach erosion at the eastern end. Furthermore, once suspended, the sediments would be transported with nearshore currents. Considering that the wave-induced currents play a dominant role in the longshore and cross-shore transport, the tidal water level applied in the simulation only varies with time throughout the model domain. However, the neglection of tidal currents still has potential impact on the beach evolution, although they only have a weak intensity.
(2) There are three sediment patterns in the model domain, and in a realistic coastal environment, these sediments could be mixed especially under storms. Actually, during the site investigation carried out in December 2012 (Hou et al., 2013), the sediment grain size D 50 along the beach were not uniform anymore after the erosion, which is coarser at the western part (0.48 mm) and finer at the eastern part (0.36 mm), showing an intense mixing of different sediment patterns. Modeling the sediment mixing processes is extremely complicated that needs to modify the sediment information at every grid point and time step in the computation, which makes the model hard to be robust. As a result, in the present model we neglected this mixing effect, which could be another reason that induced the errors.
(3) According to Cayocca (2001), it is hard to find a set of representative wave fields that can both perfectly represent the longshore and cross-shore sediment transport processes. For this specified case, the dominant process is the strong longshore sediment transport, which is proportional to H 2.5 . Nevertheless, if we focus on the cross-shore sediment transport during the storm surge, the representative wave height should then be proportional to H 6 , which brings in a dilemma in selecting the representative wave field.
(4) Located within the intertidal zone, in natural environment the hydro-sedimentological processes in the swash zone also have potential impact on the beach morphodynamics. Established on the "phase-average" frame, the present model system is more advantageous for describing the underwater processes, and it cannot simulate the sediment transport caused by wave run-up, which is more suitable for a time-dependent sub-wave model.

Profile of cross-shore flow component
The wave-induced current system is directly responsible for the sediment transport in the surfzone, and it has a complex 3D structure. Different from a 2D case, the 3D structure is not easy to present clearly in one figure. Therefore, in the following discussions we decompose the current speed into the alongshore and cross-shore components, respectively. At each of the cross-sections, the alongshore component is defined as parallel to the shoreline, and the cross-shore component is defined as perpendicular to the shoreline. Fig. 10 illustrates the vertical distribution of the wave-induced currents at WL=3.31 m during storm surge #2, in which the positive value means offshore and negative value means onshore. It must be stressed that the data shown in the figures not only contain the traditional "undertow", but also include the cross-shore component of the longshore currents caused by the irregular bed topography. Generally, the evolution of the undertow varies with wave breaking. Balanced by the wave setup gradient in the surfzone, at the shoreward side of the breaking point the undertow direction flows onshore at the surface, forming a vertical compensation flow back offshore near the seabed. Outside the surfzone, the flow direction along the vertical profile is just opposite as above. The feature indicates two phase-averaged vertical circulations around the breaking point, and Xie et al. (2017) stated that this circulation structure may be the formation mechanism of the sandbar.
Under the action of strong storm surge, the broader surfzone greatly enhances the undertow speed (maximum speed 0.5 m/s), leading the sediments to transport offshore. However, the formation of the sandbar is not obviously evolved in this beach erosion case, because the strong longshore currents along the beach rapidly carries the sediments away to the west before they could deposit to form a bar. 5.2 Profile of alongshore-flow component Fig. 11 presents the profiles of the alongshore velocity component corresponding to the above characteristic moments. The positive value means westward while the negative value means eastward. The maximum alongshore current speed occurs at the shoreward side of the breaking point, which could reach 1.5 m/s at the surface layer (Section P2) and gradually decreases both offshore and onshore. Due to the vertical turbulent mixing effect caused by breaking waves, the vertical distribution of the alongshore currents tend to be uniform especially in the middle and surface layers, and this feature has been experimentally proved by Visser (1991). At the end of the eastern breakwater, the wave diffraction effect induces a reversed circulation (see Fig. 6d), thus around section P8 the alongshore current direction points to east in the nearshore area.
Combining the above analysis, it is concluded that both actions of undertow and longshore current have significant contributions to the present beach erosion. Once the beach sediments are entrained, the undertow sweeps the sediments to the seaside, and the strong longshore current carries the sediments to the west until bypasses the western breakwater, resulting in a continuous sediment loss of the artificial beach. Fig. 12 presents the horizontal distribution of the sediment transport volume (STV) at different characteristic mo-ments, including the alongshore and cross-shore components. Note that the calculated STV values contain both of the suspended load and bed load. Generally, the alongshore sediment transport presents a smooth distribution along the beach, which is because in the model the offshore wave direction is relatively stable (NE−ENE). The maximum longshore sediment transport occurs at the middle part of the beach. In the modeled 5 months, the longshore STV caused by two storm surges takes a percentage of approximately 60%. Despite occupying most of the time, the longshore STV under calm weather only takes a percentage of 40%. During the time between two storms the longshore STV increases, the reason is that most of the coarse sands at the beach surface (D 50 =0.67 mm) have been significantly eroded after Storm #1, thus the fine sediments (D 50 = 0.13 mm) exposes, which could be easier re-suspended even under calm waves.

Sediment transport feature
The distribution of the cross-shore STV presents the feature of higher intensity in the western parts of the beach, which is because apart from the undertow, the longshore current also has a seaward deflection when bypassing the western breakwater tip (see Fig. 7) and enhances the offshore sediment transport. Under calm waves, the eastern beach (approx. 400 m) has a slight onshore sediment transport. In fact, though not significant, the calm waves would induce the net onshore drift that helps the beach to recover (Xie et al., 2017), and again the wave-induced circulation near the eastern breakwater also contributes to the onshore  Fig. 12. Distribution of longshore and cross-shore sediment transport volume.

22
XIE Ming-xiao et al. China Ocean Eng., 2021, Vol. 35, No. 1, P. 12-25 sediment transport in that region. Based on the model results, the total sediment loss volume of the artificial beach is 19.1×10 4 m 3 . Compared with that analyzed from the field measurement data (22.1×10 4 m 3 ), the model slightly underestimates the sediment loss of about 13.6%, which is in an acceptable range regarding the complexity of the specified case.
For shedding more lights on the sediment transport feature, we calculate the sediment transport rates during different characteristic periods, especially with a separation of bed load and suspended load, see Table 4. Several conclusions are drawn as follows.
(1) According to the model results, the sediment transport rate has a small magnitude under calm wave environment with only (0.03−0.04)×10 4 m 3 /d, but greatly enhances during storms with a value of (2.48−6.29)×10 4 m 3 /d, which reaches several hundreds of times as those under calm waves, showing that the short-duration storm surges contribute most of the beach erosion.
(2) From the aspect of temporal distribution, the sediment transport during Storm #2 (S2) is significantly larger than that of Storm #1 (S1). That is because not only the second storm has larger intensity, but also the exposure of fine sediments underneath (0.13 mm) are more easily suspended.
(3) Under calm waves, the sediment transport pattern is mainly bed load. During the initial period of C1, the beach is covered with coarse sands (0.67 mm), which could not be substantially suspended beneath small waves, hence the bed load transport takes a percentage of 89.2%, with suspended load transport being only 10.8%. However, as the beach erosion develops, the fine sands begin to expose. As a result, the contribution of suspended load increases to 71.1%−75.6%.

Comparison with 2DH model result
Considering that the application of 3D model has not been substantially reported, it is helpful that we compare the model results with a 2DH model. In fact, in earlier works, Zhang et al. (2019) have modelled the same case as in this paper but using a well-acknowledged XBeach model. Table 5 lists the model results of beach erosion volumes from their model and the present model. It can be visualized from the table that both models have similar results, where the 2DH model overestimates the erosion volume by 12.9%, while the 3D model underestimates the erosion volume by 13.6%, which are quite close, showing that both models could effectively reproduce the strong beach erosion, at least in this specific case.
It should be specially stressed here that although the 2DH model is commonly regarded as the depth-integration of a 3D model, there are still some essential differences between these model systems including (1) the dynamic factors considered (e.g. surface roller, net drifts), (2) the frame of sediment transport theory (e.g. reference concentration, transport rate formulation), and (3) the treatment of boundary conditions (e.g. bottom shear stress and bed roughness). Therefore, the deviation of model results between two models are not merely caused by the "dimension difference" itself, but also influenced by a series of other factors.
However, the comparison between these two model systems indicates that the 2DH model could also effectively grasp the general feature of beach evolution, which is not a surprise because for this specified case, the longshore sediment transport takes a dominant role, whose magnitude is tens of times of the cross-shore sediment transport, even under calm waves (see Fig. 12). That is because the strong longshore currents induce a continuous westward sediment Note: The time span of characteristic periods C1−C3 and S1−S2 are indicated in Table 3. flux, which cut the sediment supply from offshore, and subsequently suppress the onshore sediment deposition, as well as the formation of sandbars. Therefore, the beach "recovery" under calm waves is not obvious, which is not the case as stated by Xie et al. (2017) by analyzing a flume case of sandbar migration. That means, the beach recovery under calm waves does not necessarily exist, but with a competing process of onshore drift, undertow and longshore currents, which makes the nearshore sediment transport process extremely complex. Considering that a 2DH model could sufficiently provide the planar distribution of longshore currents, we could see in particular from Fig. 11 that the vertical distribution of longshore currents does have a near-logarithmic profile because generally the depth contour of the beach is smooth without a sharp variation, which further reduces the formation of 3D rip currents.
Generally, the 3D model is more advantageous in description of the cross-shore sediment transport due to a more sophisticate consideration of the vertical processes e.g. the undertow and sediment concentration profile, hence the good agreement between these two models (2DH and 3D) further proves the validity of the established 3D model system. If we take the model result of this paper and that in Xie et al. (2017) as an integrated validation, we could see that the 3D model could effectively reproduce a variety of beach evolution cases including the sandbar profile migration under normally incident waves (dominated by vertical undertow and drifts, cannot be modelled using a 2DH model), and also the present case (dominated by longshore currents, mainly 2DH), showing that the present 3D model is a general system, which can be used in any circumstance, especially the engineering application cases under either crossshore dominated or alongshore dominated, or cross-shorealongshore mixed.

Conclusions
A well-established three-dimensional phase-averaged beach morphodynamic model was applied to investigate the morphological performance of Weifang Artificial beach, China. Based on the model results, a series of discussions are made on surfzone hydro-sedimentological processes under calm and storm circumstances. Some key conclusions are drawn as follows.
(1) The established 3D model has successfully presented the movement of the wave propagation, wave-induced current structure as well as the beach morphology evolution. The modeled beach elevation variation agrees well against the field survey data, proving satisfactory behavior of the model for predicting nearshore hydro-sedimentological processes under realistic coastal environment.
(2) According to the model results, the nearshore waveinduced current presents a significant 3D structure under storms, where the undertow and longshore currents exist simultaneously, forming a spiral-like circulation system in the surfzone. Though the near-bed undertow swashes the sediment offshore, the strong longshore currents shorten the sediment supply in the cross-shore direction, and then suppress the formation of sandbars, showing that a typical recovery profile under calm waves does not necessarily develop once the longshore sediment transport is strong and continuous enough.
(3) Sediment transport rate during storms reaches several hundreds of times as those under calm waves. Two storm events contribute approximately 60% to the beach erosion. Sediment transport pattern under calm waves is mainly bed load, with a maximum percentage of 89.2%. As the beach erosion develops, the fine sands begin to expose, and the contribution of suspended load increases to 71.1%−75.6%.