3DModeling andMechanism Analysis of BreakingWave-Induced Seabed Scour around Monopile

Breaking wave-induced scour is recognized as one of the major causes of coastal erosion and offshore structure failure, which involves in the full 3D water-air-sand interaction, raising a great challenge for the numerical simulation. To better understand this process, a nonlinear 3D numerical model based on the open-source CFD platform OpenFOAM®was self-developed in this study. -e Navier–Stokes equations were used to compute the two-phase incompressible flow, combining with the finite volumemethod (FVM) to discretize calculation domain, a modified VOF method to track the free surface, and a k − ε model to closure the turbulence. -e nearshore sediment transport process is reproduced in view of shear stress, suspended load, and bed load, in which the terms of shear stress and suspended load were updated by introducing volume fraction. -e seabed morphology is updated based on Exner equation and implemented by dynamic mesh technique. -e mass conservative sand slide algorithm was employed to avoid the incredible vary of the bed mesh. Importantly, a two-way coupling method connecting the hydrodynamic module with the beach morphodynamic module is implemented at each computation step to ensure the fluid-sediment interaction. -e capabilities of this model were calibrated by laboratory data from some published references, and the advantages/ disadvantages, as well as proper recommendations, were introduced. Finally, nonbreakingand breaking wave-induced scour around the monopile, as well as breaking wave-induced beach evolution, were reproduced and discussed. -is study would be significantly helpful to understand and evaluate the nearshore sediment transport.


Introduction
As ocean wave gets closer to the coastal region, the free surface deforms and wave steepens with increasing wave height and may induce unstable and break waves. Wave breaking and the resulting current are the causes of many complex phenomena in the surf zone, such as the coastland drift, beach erosion, sediment transport, and so on. ey will erode the soil-supporting pile foundations [1], which were constructed within aforementioned coastal high-hazard areas for practical reasons [2]. erefore, assessing the hydrodynamic and sediment transport processes are important for inferring the strength of potential destruction during the extreme environment (e.g., Wu et al. [3,4]).
Over the last few decades, some research studies were performed to improve the understanding of above processes. Field surveys were performed by analyzing sediment diameter distribution, instrument data, and field videos (e.g., Paris et al. [5]; Szczuciński et al. [6]; FEMA [2]; and USGS [7]), which provided some substantial cognitions, such as the majority of the deposits are from beach and not from deep ocean floor (e.g., Morton et al. [1] and Jaffe et al. [8]). However, the field research studies raise considerably limited information due to frequent lack of the wave condition and topography, especially facing the complex three-dimension high-speed water roller. Additionally, theoretical approaches for studying this process are still inadequate [9]. Hence, numerical simulations and physical experiments, providing effective controlling factors, become the preferred methods.
Traditionally, physical experiments are clearly helpful to assess this process. Many research studies have been performed in the aspect of cross-shore sediment transport on sandy beach (e.g., Kobayashi and Lawrence [10]; Young et al. [11]; Jiang et al. [12]; and Daghighi et al. [13]) and the scour around offshore structure foundation (e.g., Kato et al. [14,15]; Tonkin et al. [16]; and Kuswandi et al. [17]), which improve our understanding of beach evolution and its response to structure by wave. With the development of computer technology, the numerical model provided us with another way to understand the mechanisms, which solves the limited capability used in measurement devices. In the past decades, modeling sandy beach evolution has been studied with depth-averaged equations, such as shallow water equations (e.g., Simpson and Castelltort [18] and Pritchard and Dickinson [19]) and Boussinesq-type equations (e.g., Shimozono et al. [20] and Xiao et al. [21]). However, this process is highly nonlinear, as well as local but strong turbulence near the free surface, and the seabed needs to be considered. Recently, Nakamura and Yim [22], Li et al. [23], and Jacobsen et al. ( [24][25][26]), as well as Liu et al. [27], introduced a volume of fluid-type model to study this process, which was based on the generalized Navier-Stokes equations with a turbulence closure solver computing incompressible viscous multiphase flows. Aforementioned research studies not only substantially improved our understanding of breaking wave-induced longshore sediment transport mechanisms but also gave us some inspiration to explore the more complex problems, such as sediment-wave-structure interaction in some coastal areas. Pan and Huang [28] are the pioneers to numerically investigate the interactions among wave, monopile, and sandy beach to access the effectiveness of a linked 2D hydrodynamic and sediment scour model in this process modeling. However, some limitations existed, such as the 2D model could only simulate vertically averaged currents instead of 3D currents in this study, in particular about the vertical velocity component [28]. Huang et al. [29] pointed out that the approximation of the 3D scour phenomenon around the pile structure to a 2D problem may be a contribution to errors between the model predictions and observations. To the authors' knowledge, the research of 3D modeling in this question is still limited so far and yet to be further investigated.
To address the aforementioned issues, a 3D numerical model was developed based on the open-source Open-FOAM platform in this study, supporting two-phased incompressible flow. e wave dynamic is achieved by its solver interFoam combined with the active wave generation and absorption boundary conditions, developed by Higuera et al. [30]. Unfortunately, the official version of OpenFOAM lacks the function of sediment transport simulation, and the third-party sediment module could not be open accessed. Additionally, similar sediment solving modules did not necessarily apply to current research cases due to the complexities of sediment transport mechanisms under different sea environments. erefore, this study aimed to overcome those gaps by self-extending the sediment functionality based on an effective usage of the better parts from published sediment modules.
is developed model will be further used to simulate the breaking waveinduced scour around the monopile, as well as investigate its difference with the scenario of monopile on a flat seabed and the potential scour mechanisms. e generation, advection, and dissipation of turbulence in the fluid model are also important for sediment transport. For breaking wavedriven sediment motion, several turbulence models have been employed, namely, the RANS model (Li et al. [23] and Jacobsen [25,26]) and the LES model (Christensen and Deigaard [31] and Christensen [32]). Generally, LES was under consideration as the turbulence model because of a considerable part of the mixing is resolved. However, when using the LES, the stochastic nature of the flow reflects nonlinearly onto the sediment transport and results in larger spatial gradients; hence, the allowable calculation stability or time step will be lowered considerably. us, the RANS turbulence model was considered in this study because it is more preferable to predict the scour profile than the LES turbulence model [33,34]. e rest of this paper is organized as follows. e numerical methods including the hydrodynamic model and beach morphodynamic model, as well as corresponding boundary conditions and numerical schemes are introduced in Sections 2-4, respectively. e validations of the main modules are conducted in Section 5. e model applications and discussions are given in Section 6. e main conclusions of this study are shown in Section 7.

Hydrodynamic Model
3D RANS equations were applied to simulate the two-phase incompressible viscous fluids, and the governing equations are as follows: Continuity equation: Momentum equation: where x � (x i , x j , x k ) is the Cartesian coordinate system, u � (u i , u j , u k ) is the velocity vector, ρ is the fluid density, p * is the pressure in excess of hydrostatic, and g is the gravitational acceleration. e last term of equation (2) is used to describe the surface tension, in which σ T is the surface tension coefficient, κ c is the surface curvature, and c is a scale field used to track the fluid movement. τ * is the Reynolds stress tensor given by 2 Mathematical Problems in Engineering where μ is the eddy viscosity, S � (∇u + (∇u) T )/2 is the rate of strain tensor, k is the turbulent kinetic energy, and I is the Kronecker delta function confirmed as 1 if i equal to j and as 0 if i unequal to j. e second order standard k − ε turbulent model was employed in this study to closure the set of RANS equations: where C μ is the dimensionless coefficient, l is the turbulence length scale, and ε is the turbulent dissipation rate. e transport equations of k and ε are as follows: where ] t is the turbulence kinematic viscosity and C 1ε , C 2ε , σ k , and σ ε are the empirical coefficients confirmed as 1.44, 1.92, 1, and 1.3, respectively. e free-surface motions were captured by the modified VOF method and defined as zα zt where the air and water are described by volume fraction (α) ranging from 0 (air) to 1 (water). Compared with the traditional VOF method, the extra compression term was introduced in the last term of equation (6) to limit the excessive numerical diffusion and the smearing of the interface, where u r is the relative velocity.

Beach Morphodynamic Model
In this study, the details on each component of the beach morphodynamic module are selected and tested according to a sequence of classic examples, and meanwhile, some proper parameters are recommended, which is a comprehensive usage of some previous research studies. e bed shear stress, linking the hydrodynamic features and sediment transport, is estimated by the method proposed by Arzani et al. [35]: where τ t is the wall traction confirmed as σ · n, n is the unit normal vector that is perpendicular to the bed surface mesh, and σ is the stress tensor given by where the first term is the pressure-induced stress term, p is the pressure, and the second term is the viscous stress term controlled by the bottom fluid motion. e volume fraction was introduced into equation (7) to identify the bed shear stress of the beach surface.
e bed load transport was described by the formula proposed by Engelund and Fredsøe [36]: 50 , where q b is the bed load transport rate per unit width, d 50 is the median sediment diameter, R is the relative density of the sediment, θ is the Shields number, and θ c is the critical Shields number considering the effect of seabed slope (Allen [37]) and defined as Equation (11) indicates that θ c is adjusted to a lower value for the sediment moving down the slope (i.e., negative slope) and a higher value for the sediment moving up the slope (i.e., positive angle), and φ is the sediment repose angle. θ co is threshold shields number under flat bed, which could be calculated by [38] where D * is the dimensionless sediment size.
e bed load transport rates in different directions are where η is the bed elevation and C is used to reveal the effect of bed slope on the sediment flux and is specified as 1.5 in this study. e suspended load transport can be described according to the classical convection-diffusion equation as follows: where c is the concentration of the suspended sediment, v d is the sediment diffusivity equaling to the turbulence eddy Mathematical Problems in Engineering viscosity, σ c is the turbulent Schmidt number, which is related to eddy viscosity and sediment diffusivity and is taken to be 0.8, and w s is the sediment setting velocity, which will be affected by the suspension concentration where ζ is the suspended sediment size related constant and is specified to be 5.0 in the present research according to the study of Liang et al. [33] and w s0 is the settling velocity in clear water given by [39] w s0 � v d 50 10 As far as we know, the coastal zone is the interface area of water, air, and sediment, resulting in a frequent interaction. To keep the sediment to drop out immediately when accidentally left in the air, the volume fraction was introduced in equation (15) as follows: where the velocity is multiplied by α, and the inclusion of v in equation (18) is used for numerical stability. is stability problem occurs when v d becomes 0; hence, the presence of v ensures that an advection-diffusion problem is solved locally rather than an advection problem, where the latter is more difficult to solve numerically [39]. According to sediment continuity, the Exner equation was used to describe the seabed evolution: where z is the seabed elevation, n is the beach porosity, and q b is the bed load transport rate, and its components are given by equation (14). D is deposition rate solved by where c s is the sediment concentration at bed defined as the value at the nearest cell center of beach in this study. Additionally, the sediment concentration at a reference level Δ b is calculated by the following formula: where T is the dimensionless excess shear stress calculated by in which τ cr is the critical shear stress and μ sc is the effective shear stress coefficient. e details could be found in van Rijn [40]. E is the entrainment rate, which could be calculated by the following formula: Finally, the model coupling between water and sediment is achieved by moving the computational mesh in such a way that the bottom mesh is conformal with the seabed, which is mainly based on the method proposed by Jasak and Tukovic [41]. Meanwhile, to ensure that the bed slope does not exceed the sediment repose angle, the mass conservative sand slide algorithm given by Khosronejad et al. [42] was employed: where z bp and z bi are the bed elevations at point P and its ith neighbor (see Figure 1), Δl pi is the horizontal distance between the two cell centers, and Δz bp and Δz bi are the corresponding corrections imposed to satisfy the sediment repose angle, which could be obtained by balancing the cell mass as follows: where A bp and A bi are the projections of the cell P and its ith neighbor.

Boundary Conditions. Proper boundary conditions are
the key factors to copy laboratory experiments. In this study, the active wave generation and the absorption boundary developed by Higuera et al. [30] were employed. Compared with the traditional method proposed by Jacobsen et al. [24], this method does not need the relaxation zone so that it could significantly reduce the computational domain. e other boundary conditions are set as follows. e top boundary is described as free to the atmosphere, and the seabed and the pile surface are set as the no-slip condition. e symmetric condition is applied at the two-side faces of the numerical tank to reduce the computational cost, in view of the symmetric distribution in the streamwise direction.

Numerical Schemes.
In this study, the finite volume method (FVM) and Euler scheme are used to address the space and time term, respectively. e PIMPLE algorithm is employed for the pressure-velocity solver, which is a mixture between PISO (pressure implicit with the splitting of operators) and the SIMPLE (semi-implicit method for pressure-linked equations) method.
e MULES (multidimensional universal limiter for explicit solution) scheme is applied to maintain the boundedness of volume fraction.
e Gauss linear corrected scheme and Gauss linear scheme are used to solve the Laplacian term and gradient term, respectively. e sediment module is addressed as follows.
e bed load is solved firstly by shear stress that is derived from the hydrodynamic model. Meanwhile, the suspended load is transported as the fluid motion based on the convection-diffusion equation, and this equation is also solved by FVM. Subsequently, the two-dimensional Exner equation is applied to update bed elevation, along with the finite area method (FAM) addressing the problem of fluid information mapped from the three-dimensional space to the two-dimensional space (Tuković [43]). e FAM is a variant of the FVM, operating on two-dimensional curved surfaces in the three-dimensional space, which makes the seabed morphology suitable for arbitrary curved surfaces. In the layout, the FAM follows the structure of the Finite Volume discretization library, sharing the basic field algebra, linear system support, boundary condition settings, and discretization techniques. Subsequently, the sand slide procedure is implemented to modify the problem of bed slopes exceeding the sediment repose angle.
e self-developed model in this study is a two-way coupling scheme between water motion and sediment transport, and the computational procedure is implemented in Figure 2. Firstly, the hydrodynamic parameters, such as u n− 1 , p n− 1 , v n−1 t , α n− 1 , are computed using the hydrodynamic module. Second, above parameters are used to drive suspended load and bed load transport, resulting in the seabed change based on the Exner equation and the sand slide equation at the same (n − 1)th time step, i.e., z n− 1 . ird, the hydrodynamic parameters are also affected by the resulting seabed change. Fourthly, the values of hydrodynamic at the next nth time step are calculated by the hydrodynamic module from the previous results based on the updated bed elevation. Fifthly, the sediment transport and seabed change at the nth time step are performed using not only the previous sediment environment at (n − 1)th time step but also the hydrodynamic results at the nth time step. Finally, this process is repeated until the calculated time reaches the preassigned time. During the computation, the time step was automatically adjusted to ensure the Courant number c r (c r � Δt × max(|Δu|)/min(|Δx|), in which max(|Δu|) and min(|Δx|) were the maximum velocity and the minimum grid size, respectively, and Δt is the time step) is always less than one.

Wave Module.
Reliable hydrodynamic is the foundation to ensure the accuracy of sediment results. erefore, the stability of wave propagation over a flat tank was first tested. e expression of the solitary wave proposed by Lee et al. [44] was used to reproduce the free surface: where η is the free surface elevation, H is the wave height, h is the water depth, X � x − ct, and c � �������� g(h + H) is the wave celerity. Based on one typical solitary wave condition (H/h � 0.3), seven wave gauges with a constant interval 0.48 m were placed along the tank to obtain the wave surface elevation, and the results are shown in Figure 3. e good agreement between numerical results and theoretical results was found. In addition, the wave surface was almost no attenuation in the numerical domain, indicating the numerical model could better model the wave dispersion and nonlinearity and could be employed in the following research.

Sand Slide Module.
In this section, a constant beach slope of 45°with sediment repose angle of 30°was used to validate the robustness of the sand slide module. e beach morphology before and after employing the sand slide model is shown in Figure 4, respectively. We could find that the beach slope is less than or equal to the sediment repose angle, which indicates the sand slide module can work well once the beach slopes exceed the sediment repose angle.

Suspended Sediment Module.
e transport capacity of the suspended sediment was tested by the net entrainment experiment of van Rijn [40]. e sand bed was laid after a rigid bed in this experiment, in which x � 0 is the starting position of the sand bed, and the inlet of the model is free of sediment (see Figure 5). e flow is stable with a water depth of 0.25 m and a velocity of 0.67 m/s, and the representative sediment diameter is 0.2 mm with the corresponding settling velocity of 0.022 m/s. e roughness height k s is taken to be 0.01 m, suggested in accordance with the research of van Rijn [45] and Wu et al. [46]. e concentration profiles at x � 4H, 10H, 20H, and 40H were chosen to validate the simulated results. In this model, the entire computational domain was discretized using a structured mesh, and the gird sizes of 0.002 m were kept constant in the x-direction, y-direction, and z-direction. Figure 6 shows that the calculated results match well with the experimental data, and the sediment concentration increases with the decreasing height above the bed surface.

Bed Morphodynamic Module.
e experiment of Mao [47] was employed to validate the capacity of sediment transport with the morphodynamic module. e submarine pipeline with diameter of D � 0.1 m was laid on the sand bed, water depth is h � 0.35 m, sediment diameter is d 50 � 0.00036 m, sediment density is ρ sed � 2650 kg/m 3 , flow velocity is u � 0.35 m/s, and critical shields parameter is θ c � 0.048. To avoid the divergence of numerical results at the junction of the seabed and the pipeline caused by the severe deformation of grids, a sinusoidal hole with 0.1 D was used under the pipeline inspired by the research studies, such as Liang et al. [33], Brørs [48], and David et al. [49]. e mesh resolutions around pipelines were increased to solve the problem of mesh distortion as the scour hole broadens and deepens, as shown in Figure 7.
e scour profiles at 10 min and 30 min were chosen to verify the seabed scour morphologies, as shown in Figure 8. Overall, the good agreement between the simulation and experiment was found, in which the sediment accretion in the former behind the pipe was slightly overpredicted than that in the latter. is is mainly attributed to the fact that the RANS model tends to smooth out the fluctuations produced by the vortex shedding and therefore underestimates the interaction between the vortices shed with the seabed, and the similar conclusion was raised by Liang et al. [33].

Mathematical Problems in Engineering
Nevertheless, the sediment module established in this study could well reproduce the morphological change of the seabed. e water-sediment transport capacity of the model had been verified well by the above calculation settings. Furthermore, this model would be used to show its robustness of copying the scouring process around a pile. e numerical setup for the model validation is the same as the experiment by Sumer et al. [50]. e flume is 28.0 m long, 2.0 m wide, and 1.0 m high, and the pile diameter D � 0.10 m fixed in the center of the sand bed, as shown in Figure 9. e sediment diameter is d 50 � 0.00018 m, sediment density is ρ sed � 2700 kg/m 3 , critical shields parameter is θ c � 0.047, and water depth is h � 0.40 m. A periodic wave with wave period of T � 4.5 s and wave height of H � 0.12 m was generated to propagate over the flat sand bed. Figure 9(c) shows the model grids, and the grids around the pile were refined. Figure 10 shows the computed and theoretical wave surface elevation proposed by Svendsen et al. [51], as well as the computed and measured dimensionless maximum scour depth (S/D) change around the pile, in which S was presented without distinguishing between specific locations such as behind or beside the pile. e simulated wave crests and troughs are almost equal and in of the phase compared to the theoretical data. Meanwhile, it could be seen that the scour depth developed with the periodic wave action is well captured. Nevertheless, a small difference between the simulated and the experimental scour topography could be found. e measured scour depth fluctuated frequently and even exhibited an obvious disparity at two adjacent moments, but the variation of simulated results is rather regular, which should be attributed to the turbulence model based on a RANS method used in this study that tends to smooth out the fluctuations produced by the vortex shedding. Figure 11 shows the scouring process, and it is seen that the seabed at the upstream side of the pile was eroded after the wave impacting the pile, which is related to the strengthened circling flow. en, the suspended sediment was deposited downstream. As the wave continue impacting the pile, the scour depth at the upstream side continuously increased and the sedimentation at the downstream side also  Mathematical Problems in Engineering exhibited the same tendency. Finally, the simulated maximum scour depth took place at the upstream side of the pile, which follows the popular rule. erefore, the current model is capable of modeling the interaction between the periodic wave and the pile.

Model Setup.
e water-sediment transport capacity of the model had been verified well by the above calculation settings. In this section, the numerical evolution of sandy beach was calibrated by the experimental data conducted by Kobayashi and Lawrence [10]. e experiment setup is shown in Figures 12(a) and 12(b), with the beach slope of 1 : 12, sediment diameter of 0.18 mm, porosity of 0.4, and water depth of 0.8 m. In this experiment, a solitary wave with wave height of 0.216 m was generated to propagate over sandy beach. Another solitary wave with the identical height was generated on the evolved beach after this wave subsided, and that cycle repeats. Finally, the experimental results of beach profile, wave surface, and velocity at the end of the fourth wave were selected to test the robustness of the model. Figure 12(c) shows the model grids for sandy beach, and the grid near the seabed was refined. Figure 13 shows the comparison of wave surface elevations from the numerical model and physical experiment. Good agreements between them could be found, not only the free-surface elevations but also the arrival times. Figures 14(a) and 14(b) show the evolved beach profile and profile change. We could found that the numerical results in the offshore zone agree well with the measured data. However, the beach profile near the coastline could not be predicted accurately, which may be due to an overestimation of the turbulent dissipation because of the underestimation of the backwash velocity. Meanwhile, the time series of streamwise velocity (u) at the ADV location is shown in Figure 14(c). e numerical model satisfactorily reproduced the streamwise velocity, despite there are still some issues to be solved in underestimating the backwash velocity and the associated beach profile evolution. After completing the model validation, the wave field, suspended load transport, and the associated beach profile are shown and discussed in Figure 15. Figures 15(a)-15(c) show that the ocean wave got steepened and skewed with close to the coastal zone, and then the wave broke as a plunging breaker due to wave shoaling. e breaking waves inject nearshore immediately forming a complex wave field entraining air on the free surface, picking up amounts of   Mathematical Problems in Engineering sediment from sandy beach into suspension. Subsequently, the suspended sediment further runup, along with the breaking wave-induced current. e beach profile have not obviously changed in this stage, which should be attributed to the upwash direction opposite with the sediment gravity and bottom friction, which can also be found in the laboratory by Young et al. [11]. After upwash reaching the maximum height along the beach, it will start to reverse, as shown in Figures 15(d)-15(f ). A hydraulic jump is formed when the retreating water    Mathematical Problems in Engineering tongue collided with the relatively still massive body of water, but there is no obvious sediment resuspension. As upstream retreating water continues to inject into the hydraulic jump zone, the still water body turns to a complex water-air mixing zone that led to a turbulent roll, forming a recirculation region. Meanwhile, an intensive sediment transport in the form of sheet flow was found on the landside of the hydraulic jump point, resulting in net erosion of the onshore region. In contrast, a mass of sediments were deposited on the seaside of the hydraulic jump point due to the sudden deceleration of sediment-rich seaward flow, as well as the long particle residence time induced by the recirculation region. Consequently, the bed morphology of net erosion in the nearshore zone and net deposition in the offshore zone were formed, which was consistent with the measurements by Kobayashi and Lawrence [10] and Young et al. [11].  sandy beach, given by Tonkin et al. [16], was applied to validate the model robustness. In this experiment, the water depth was 2.45 m, wave height was 0.22 m, beach slope was 1 : 20, sediment diameter was 0.35 mm, and sediment density was 2643 kg/m 3 . A single pile with the diameter of 0.5 m was placed upright on the location of still water level as shown in Figures 16(a) and 16(b). e wave gauges were applied to measure water surface elevation at the locations of stations A, C, and D. An electromagnetic flowmeter, placing at station B, was employed to measure the streamwise velocity.

Verification and Discussion.
e scour depths at the front, side, and back of the pile, i.e., stations C, D, and E, were measured using digital camera technology. More details could be found in Tonkin et al. [16]. e numerical model was setup in the same dimension with this experiment, and the mesh is shown in Figure 16(c) with the refined grids around the pile. Figure 17 shows the time series of free-surface elevations at stations A, C, and D and the time series of streamwise velocity at station B from numerical results and measured data. Overall, the predicted results were in good agreement with the measured data, with the exception of the velocity since t � 13 s after the initial wave impact. Numerical and measured differences of Figure  17(b) after t � 13 s should be attributed to during the most turbulent part of the drawdown, at which the water level is too low to give a valid velocity measurement; the sensor is not fully submerged at this point (Tonkin et al. [16]). Fortunately, the current model is better to make up this deficiency and shows the subsequent velocity process until about t � 20 s. Nevertheless, due to the further reduction of water level, the time series of velocity start to exhibit the intense numerical oscillation.  Figure 18 shows the time series of the local scour depth around a single pile. We could find that the seabed morphology at the front and side of the pile (i.e., stations C and E) was sharply eroded of approximately 10 cm when the moment of wave impacting the pile, which is related to the stronger velocity (see Figure 17(b)). en, the scour depth decreased as the net deposition occurred during the drawdown. e seabed morphology behind the pile (i.e., station D) did not exhibit an obvious change at the moment of the wave impacting the pile, and then the scour depth continuously increased from t � 12 s to t � 18 s. During the latter stage of the backwash, the scour depth decreased to some extent due to the sediment settlement. After about t � 20 s, the seabed around the pile was relatively in equilibrium. Finally, we noted that the predicted maximum scour depth is slight lesser than the     measured data, which may be caused by seepage flow (or groundwater flow) as described in the Tonkin et al. [16], but the current model was incapable of reproducing it.

Verification and Discussion.
To more visually present the topography evolution of sandy beach during the interaction between the wave and the monopile, the three-dimensional scour patterns are shown in Figure 19. It is obvious that the remarkable erosion took place at the front of the pile when the wave impacting the pile, which is different with the aforementioned scenario of the pile placed on a horizontal seabed. As upstream retreating water, the scour depth on the seaside of the pile increased quickly. e vorticity fields, which are related to the sediment transport (Kobayashi and Oda [52]), were chosen to reveal the potential causes of above phenomena, as shown in Figure 20. At the runup stages, the flow separated and the vortex motion formed and developed mainly located at the area less than 0.5D behind the pile. Following the wave drawdown, the wake region appeared before the pile. e generated vortex was transported further upstream by advection, and the wake region was stretched with a length scale more than 2D. is may be related to the backwash velocity larger than the uprush velocity, which is showed in Section 6.2. Overall, it is acceptable that backwashinduced scour around a pile is larger than the uprush.

Conclusions
A nonlinear three-dimensional coupled model, which was composed of the hydrodynamic module and the self-developed beach morphodynamic module, was developed based on the CFD tool OpenFOAM in this study. To better simulate the breaking wave-induced seabed scour around the monopile, an innovative combination of sediment components was proposed which is different with published papers, and its robustness was verified by a series of tests. e results showed that this self-developed model was able to satisfactorily reproduce the wave transmission, suspended load and bed load motion, and the seabed morphology evolution, as well as proved the sand slide module works well. e 3D model developed here was first applied to explore two preliminary scenarios, one with the nonbreaking waveinduced scour around the monopile, and the other with breaking wave-induced sandy beach evolution. In both cases, the calculated results were basically in agreement with the measurements and meanwhile showed us that the obvious scour occurred in the upstream side of the pile under the nonbreaking wave, as well as the backrush-induced beach deformation is larger than the uprush. Subsequently, the more complicated case of wave-induced scour around the monopile on sandy beach was simulated, and its main features were well reproduced by this model. e remarkable erosion took place at the front of the pile when the wave impacting on the pile during the wave runup, which is different with the scenario of the monopile placed on a horizontal seabed. Following the wave drawdown, scour depth at the upstream side of the pile increased quickly, and this should be attributed to the high-speed backwash and strengthen the wake region. Finally, we remark that further research studies are needed to investigate the effect of the offshore structures on beach morphological evolution based on the current study, which would be significantly helpful to evaluate the nearshore sediment transport.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.