# Numerical Study of Wave Dynamics over Island Reefs Based on Smoothed Particle Hydrodynamics Method

Authors: Degang Wang; Siheng Li; Li Sun; Zhenke Zhang; Panyang Huang
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Coastal areas are suffering from an increasing number of hazards like storm surge or flooding. Luckily, reefs can reduce the risk because of their particular wave dynamic characteristics. In this study, a numerical wave flume is established to simulate the topography of fringing reefs around islands. The smoothed particle hydrodynamics- (SPH-) finite particle method (FPM) scheme is used to deal with Navier–Stokes equations, and the boundary repulsion method is used to deal with the solid wall boundary. The wave propagating and breaking processes over island reefs under regular wave incidence are simulated. The results are compared with the experimental results in the literature, showing that the numerical model can simulate the wave propagating process on the island reef terrain very well. The wave propagation over the island reef terrain experienced the processes of shoaling, breaking, reforming, and so on, and most waves break in the plunging form. The processes of wave-induced setup accompanied with wave breaking before and after the reef ridge were also simulated.


---

## Body

## 1. Introduction

Large numbers of people live near the coast. Unfortunately, researchers found that wind speed, wave height, and extreme events are increasing globally in the sea [1], threatening human life. Reefs have proved to be effective in coastal hazard risk reduction. Different landforms from the traditional coast will change wave dynamics. 97% of the wave energy could be dissipated by reefs [2], avoiding the waves to impact the coastline. In addition, reefs have other effects because of their special topography. For example, in the nearby waters, the nearshore circulation promotes the movement of nutrients, which are crucial to the production, dispersal, and retention of larval fish and other invertebrates [3]. Meanwhile, sediment transport depends on water movement associated with waves in the reef area, which affects shoreline stability a lot [4].In the 1950s, researchers began a series of field observations to investigate wave propagation characteristics. Munk and Sargent [5] noticed that mean water levels in the surrounding ocean tend to be lower than those on the reef flat of Bikini Atoll. Based on data from field observation, Tait [6] found that the wave-induced setup can be predicted well if momentum flux is considered. When a wave transports over the reef, at least 85% of its energy will be absorbed at both high tide and low tide [7]. Excitation of the resonant modes can be aroused due to the wave setup. Subsequently, resonant frequencies will increase and result in the occurrence of wave group forcing with time scales as a result [8]. Meanwhile, many laboratory studies focused on wave breaking and wave-induced setup over reefs as well. Wave breaks on the reef crest and the reef-face or the reef-rim at high tide and low tide, respectively. Relevant theories are radiation stress theory and the continuity theory, which can predict corresponding wave setup well [9]. There exist relationships between wave setup and wave-generated flow, which can be utilized to predict experimental values in most cases [10]. Based on their experimental data, Yao et al. put forward empirical formulas in order to investigate wave runup [11]. In addition, they found that permeability of reef flat is a significant factor in the wave breaking and attenuation by comparing transportation of regular waves over smooth impermeable and rough permeable reef models [12].Theoretical models and numerical studies are important supplements to physical model tests. Schaffer et al. [13] developed a model, which can describe the occurrence and stop of wave breaking and the evolution of breaking wave profile. By assuming friction laws and wave breaking criteria, wave-generated current over the reef was investigated as well as its relationship with different topography and roughness of the reef [14]. Based on the mild slope equation and the radiation stress theory, the intensity of wave setup and wave-induced current were estimated by Massel and Gourlay [15]. Yao numerically investigated the tsunami-like solitary waves propagating over certain photography [16]. Similarly, the Boussinesq model has proved to be capable of predicting the wave height and mean water level over different reef profiles and describing the complex changes to the wave spectrum over reef flat [17, 18].Mesh-based methods, including Finite Volume Method (FVM), Finite Difference Method (FDM), and Finite Element Method (FEM), have been widely applied in numerical simulation. However, there exist some problems that the mesh-based method is hard to handle, such as problems involving complex topological changes, large deformation, or moving interface [19]. Smoothed particle hydrodynamics (SPH) has been developed in recent years, which is suitable to deal with the abovementioned challenges. It is a meshless particle method with strong self-adaptability [20]. Kernel approximation and particle approximation are two main steps in the SPH scheme, after which field variables on a particle are approximated with a weighted average of the values over neighboring particles. However, the conventional SPH method still needs much improvement, such as its particle consistency and uniform distribution. For example, the δ-SPH scheme was put forward by Autuono et al. [21, 22] to overcome spurious pressure fluctuations. In addition, they proposed an arbitrary Lagrangian-Eulerian framework based on the standard differential formulations of the WCSPH [23]. Xu [24] developed particle shifting technology to achieve better particle distribution. Huang et al. [25] proposed the coupled SPH-FDM in an Euler frame for unsteady flows around the pitching airfoil. They not only expanded the application of the SPH method but also gave some significant conclusions for the study of pitching airfoil. To accelerate simulation and overcome a large quantity of calculations in problems with large scales, Wang et al. [26] utilized a graphics processing unit (GPU) instead of the central processing unit (CPU) to obtain better parallelism. The finite particle method (FPM) proposed by Liu et al. [27, 28] is a modified SPH scheme with better accuracy [28]. In recent years, the SPH and FPM scheme have been applied to solving different coastal engineering problems, such as wave breaking and overtopping [29–35], fluid-structure interactions [36–40], and porous media [41–47].This study establishes a 2D numerical model to simulate how regular waves transform, break, and lead to wave-induced setup over a fringing reef and how the reef ridge and wave parameters influence this process. FPM scheme coupled with SPH is applied to solve the Navier–Stokes equations, which play the role of governing equations. As to the boundary condition, a modified boundary force model with a repulsive force on the fluid particles adjacent to the boundary is used to deal with solid wall, while dispersion relation and wavemaker theory are applied additionally to the piston-type wavemaker. The numerical results are in good agreement with experimental results obtained by Yao et al. [48], which proves that the SPH-FPM scheme can simulate the wave transformation process over reefs accurately. There are four sections in this paper. In Section 2, the methodology is introduced, including governing equations, coupled SPH-FPM scheme, and boundary conditions. In Section 3, the computational model is described, and then the regular waves, the wave breaking process, and wave-induced setup are discussed, respectively.


## 2. Methodology

### 2.1. Governing Equations

For a weakly compressible media, the continuity and momentum equations in a Lagrangian framework are shown as follows:(1)dρdt=−ρ∇⋅u,(2)ρdudt=−∇p+τ+F.In equations (1) and (2), ρ is the initial density, u is the velocity vector of a particular position, p is pressure, τ is viscous force, and t is time. F is the external body force, which is equal to the gravity g and thus set to be 9.81 m/s2. ∇ and ∇ are the divergence and the Hamiltonian operators, respectively. The pressure is calculated by the following equation (3) in the weakly compressible state [49]:(3)p=c2ρ−ρ0,where c is the artificial sound speed instead of actual sound speed to avoid an unacceptably small time step. ρ0 is the pressure of free surface. Here, the artificial sound speed is dictated by(4)c=10gh,where g is the gravity force and h is the water depth before the fore-reef slope in the beginning.


### 2.2. SPH Scheme

In the SPH method, the whole space is represented by a finite number of particles with independent mass and independent space occupied. There are two critical steps in the SPH method. Firstly, in the process of kernel approximation, the kernel functionW is introduced to denote the interaction among SPH particles. Secondly, for a given parameter of a particle, the kernel function and its derivative are replaced by the sum of the corresponding values of adjacent particles in its support domain, which is called particle approximation. By applying kernel approximation and particle approximation, the approximate formulas of a function and its derivatives are as follows:(5)fxi=∑j=1NmjρjfxjWij,(6)∇⋅fxi=∑j=1Nmjρjfxj⋅∇iWij,(7)Wij=Wxi−xj,hs=Wxi−xj,hs,(8)∇iWij=xi−xjrij∂Wij∂rij=xijrij∂Wij∂rij,where  represents the SPH approximation, x is the spatial coordinates, xj is the spatial coordinates of interacting particles within the support domain of xi, and h is the smoothing length determining the support domain of kernel function. W (xi−xj, hs) is the kernel function for the relative position xi−xj and smoothing length hs. N is the number of particles within the support domain of particle i, and rij is the distance between particle i and particle j. The cubic spline kernel function is adopted as(9)WR,hs=αd×23−R2+12R3,0≤R<1,162−R3,1≤R<2,0,R≥2,where αd is a coefficient with the value of αd=15/7πhs2 in a 2D space and R = rij/hs.With equations (5)–(8) applied to equations (1) and (2), the continuity and momentum equations can be expressed as [20, 50](10)Dρidt=ρi∑j=1Nmjρjuij•∂Wij∂xi,ρiduidt=−∑j=1Npj+pi∇iWVj+τi+Fi,where Vj is the volume of the particle, Vj=mj/ρj, and uij=ui−uj. The following artificial viscous term t can be used in the context of viscous free-surface flows [51]:(11)τi=αhscρ0∑juj−ui⋅rj−rirj−ri2∇iWVj,where α is the artificial viscous coefficient, which proves to be of great significance to the simulation and usually in the range [0.01, 0.05] [52]. And, α is set to be 0.05 in this study.


### 2.3. Improved SPH Scheme with the FPM Method

However, the conventional SPH method is limited in approximation accuracy despite its favourable stability, probably because of particle inconsistency. As an implementation of the conventional SPH scheme, Liu et al. [27, 28] proposed a finite particle method (FPM), which restores particle consistency and thus improves the approximation accuracy. By applying standard FPM, the continuity and momentum equations can be written as(12)dρidt=−ρia21∑jujWVj+a22∑juj∇i,xWVj+a23∑juj∇i,yWVj,(13)ρiduidti=−a21∑jpjWVj+a22∑jpj∇i,xWVj+a23∑jpj∇i,yWVj+τi,x+Fi,x,ρidvidti=−a31∑jpjWVj+a32∑jpj∇i,xWVj+a33∑jpj∇i,yWVj+τi,y+Fi,y,where u and v are tangential and normal velocities, respectively, and ∇i,xWand∇i,yW are kernel derivatives with respect to the particle i in x and y directions. Coefficients ααβα,β=1,2,3 are components of corrective matrix A [53]:(14)A=a11,a12,a13a21,a22,a23a31,a32,a33=∑jWVj∑jxjiWVj∑jyjiWVj∑j∇i,xWVj∑jxji∇i,xWVj∑jyji∇i,xWVj∑j∇i,yWVj∑jxji∇i,yWVj∑jyji∇i,yWVj−1.Given a well-conditioned matrixA, FPM will have higher stability and better accuracy. So, we want particle distribution to be as uniform as possible during the simulation process. Particle shifting technology (PST) can shift particles slightly through streamlines every step, thus alleviating particle clustering or sparse discretization. The PST includes the following two steps [24]:Step 1: use the following formula to modify the particle position:(15)δri=−hs22∑j1+0.2WWΔ4∇iWVj,whereδri is the displacement correction and WΔ is the value of the smoothing function in the position of the particle itself.Step 2: after completing step 1, use the following formula to correct the field variables of particles:(16)φi′=φi+δri⋅∇φi,wherei′ and i represent modified field variables and previous ones and φ stands for a general field variable.As expressed above, the FPM scheme has more advantages compared to SPH. However, FPM is usually not suitable for particles near the free surface because condition numbers here tend to be very large. At the same time, the SPH scheme naturally satisfies the dynamic boundary condition of the free-surface flows. Therefore, in this study, the SPH scheme is still applied to calculate the field variables of fluid particles near the free surface, whereas the FPM scheme is applied to the rest. Particles in the vicinity of the free surface are identified as free-surface particles and particles neighboring them, as shown in Section3.2.Here we talk about the criterion of whether a particle is located near the free surface at a certain time step [54]. First, we take the particle as the origin, make the horizontal and vertical coordinate axes, and divide its support domain into four parts (i.e., the radius of the blue circle in Figure 1 is khs), in which the number of fluid particles is Q1–Q4, respectively. Then, the coordinate axis is rotated π/4, with the number of particles in the four quadrants recorded as Q5–Q8. For the fluid particles in the vicinity of the free surface, at least one of Q1–Q8 is less than the constant K1. Otherwise, it is an interior particle. B-spline kernel function is used as the kernel function, and the corresponding constant K1 is equal to 3.Figure 1
Illustration of the method for detecting free-surface particles.
(a)(b)


### 2.4. Boundary Conditions

The particles include fluid particles and virtual particles. Virtual particles are the physical boundary in the simulation, preventing the fluid particles from passing through the boundary and providing necessary information for fluid particles in its support domain. The virtual particles involved in the simulation process can be divided into two categories according to their positions. One is virtual particles to represent the piston-type wavemaker, and the other one represents fixed solid boundary. The horizontal velocity and displacement of the virtual particles representing the wavemaker are calculated by the wave-making theory, while the position of the virtual particle forming the fixed boundary remains stable, and its velocity is determined by the velocity at the projecting point.


#### 2.4.1. Solid Wall

In the following simulation, the projecting points associated with the boundary particles are created to calculate the physical information of boundary particles as shown in Figure2. The projecting point xp is determined by(17)xp=xb+εnbΔ,where xb is the position of the virtual particle relevant to the projecting point [55], nb is the normal vector of the boundary particle, and Δ is the initial particle spacing. ε is the factor to make sure that the projecting point locates at half the initial particle spacing away from the boundary line.Figure 2
Diagram of projecting points corresponding to boundary particles.The fluid variables (velocity, pressure, and density) of the projecting point can be evaluated by(18)fp=∑mj/ρjfjWpj∑mj/ρjWpj,where j represents all fluid particles in the support domain of the projecting point.Field variables of virtual particles are calculated differently from those of fluid particles. Velocity vectors are maintained by(19)ubτ=upτ,ubn=−upn,where ubτ and upτ are the tangential velocities of virtual particles and projecting points. un represents the normal velocity. Pressure and density should satisfy the weakly compressible equation as(20)pb=pp+ρ0ghb−hp,(21)ρb=ρ0+pbc2,where pb and ρb are the pressure and density of virtual particles, ppandρp are the pressure and density of projecting points, and hbandhp are the water depth of virtual particles and projecting points.In the process of numerical study, fluid particles may get close to or even pass through the boundary formed by virtual particles. To prevent this kind of error, virtual particles should exert a repulsive force on fluid particles in their support domain as(22)fij=−cui−uj∗njWijhij2njrij∗nj,ui−uj∗nj<0,0,ui−uj∗nj>0,where rij=xi−xj and hij=hsi+hsj/2, with nj the normal vector and c artificial sound speed.


#### 2.4.2. Wavemaker Theory

Wave parameters are not independent of each other and should satisfy the dispersion relationship. Given water depthh, wave height H0 and wave period T are taken as variables in the following simulation process, and other wave parameters can be determined by the dispersion relationship as(23)σ2=gKtanhkh,where σ is the angular frequency, which is dictated by(24)σ=2πT.Thus, phase speedc0 and wavelength L are then available:(25)c0=gT2πtan hkh,L=c0T.The incident wave in our simulation is generated by a piston-type wavemaker located on the left side of the flume. The horizontal displacement and velocity of the piston are given, respectively [56]:(26)xpistont=s2sinσt,upistont=σs2cosσt,where s is defined by s=H0/kh. Here, H0 is the wave height, h is the water depth, and k is the wavenumber.


### 2.5. Time Integral Scheme

The Predictor-Corrector time integration algorithm is used here. At the first half step, all fluid variables are predicted using the following schemes:(27)ρin+1/2=ρin+Δt2dρindt,vin+1/2=vin+Δt2dvindt,xin+1/2=xin+Δt2vin.Then, these values at the second half step can be calculated as(28)ρ¯in+1/2=ρin+Δt2dρin+1/2dt,v¯in+1/2=vin+Δt2dvin+1/2dt,x¯in+1/2=xin+Δt2vin.Finally, all fluid variables at the next time step can be calculated as follows:(29)ρin+1=2ρ¯in+1/2−ρin,vin+1=2v¯in+1/2−vin,xin+1=2x¯in+1/2−xin.To ensure the accuracy and stability, the time step is selected by(30)dt=CFL∗Δc.Herein, the CFL number is chosen as 0.5.


## 2.1. Governing Equations

For a weakly compressible media, the continuity and momentum equations in a Lagrangian framework are shown as follows:(1)dρdt=−ρ∇⋅u,(2)ρdudt=−∇p+τ+F.In equations (1) and (2), ρ is the initial density, u is the velocity vector of a particular position, p is pressure, τ is viscous force, and t is time. F is the external body force, which is equal to the gravity g and thus set to be 9.81 m/s2. ∇ and ∇ are the divergence and the Hamiltonian operators, respectively. The pressure is calculated by the following equation (3) in the weakly compressible state [49]:(3)p=c2ρ−ρ0,where c is the artificial sound speed instead of actual sound speed to avoid an unacceptably small time step. ρ0 is the pressure of free surface. Here, the artificial sound speed is dictated by(4)c=10gh,where g is the gravity force and h is the water depth before the fore-reef slope in the beginning.


## 2.2. SPH Scheme

In the SPH method, the whole space is represented by a finite number of particles with independent mass and independent space occupied. There are two critical steps in the SPH method. Firstly, in the process of kernel approximation, the kernel functionW is introduced to denote the interaction among SPH particles. Secondly, for a given parameter of a particle, the kernel function and its derivative are replaced by the sum of the corresponding values of adjacent particles in its support domain, which is called particle approximation. By applying kernel approximation and particle approximation, the approximate formulas of a function and its derivatives are as follows:(5)fxi=∑j=1NmjρjfxjWij,(6)∇⋅fxi=∑j=1Nmjρjfxj⋅∇iWij,(7)Wij=Wxi−xj,hs=Wxi−xj,hs,(8)∇iWij=xi−xjrij∂Wij∂rij=xijrij∂Wij∂rij,where  represents the SPH approximation, x is the spatial coordinates, xj is the spatial coordinates of interacting particles within the support domain of xi, and h is the smoothing length determining the support domain of kernel function. W (xi−xj, hs) is the kernel function for the relative position xi−xj and smoothing length hs. N is the number of particles within the support domain of particle i, and rij is the distance between particle i and particle j. The cubic spline kernel function is adopted as(9)WR,hs=αd×23−R2+12R3,0≤R<1,162−R3,1≤R<2,0,R≥2,where αd is a coefficient with the value of αd=15/7πhs2 in a 2D space and R = rij/hs.With equations (5)–(8) applied to equations (1) and (2), the continuity and momentum equations can be expressed as [20, 50](10)Dρidt=ρi∑j=1Nmjρjuij•∂Wij∂xi,ρiduidt=−∑j=1Npj+pi∇iWVj+τi+Fi,where Vj is the volume of the particle, Vj=mj/ρj, and uij=ui−uj. The following artificial viscous term t can be used in the context of viscous free-surface flows [51]:(11)τi=αhscρ0∑juj−ui⋅rj−rirj−ri2∇iWVj,where α is the artificial viscous coefficient, which proves to be of great significance to the simulation and usually in the range [0.01, 0.05] [52]. And, α is set to be 0.05 in this study.


## 2.3. Improved SPH Scheme with the FPM Method

However, the conventional SPH method is limited in approximation accuracy despite its favourable stability, probably because of particle inconsistency. As an implementation of the conventional SPH scheme, Liu et al. [27, 28] proposed a finite particle method (FPM), which restores particle consistency and thus improves the approximation accuracy. By applying standard FPM, the continuity and momentum equations can be written as(12)dρidt=−ρia21∑jujWVj+a22∑juj∇i,xWVj+a23∑juj∇i,yWVj,(13)ρiduidti=−a21∑jpjWVj+a22∑jpj∇i,xWVj+a23∑jpj∇i,yWVj+τi,x+Fi,x,ρidvidti=−a31∑jpjWVj+a32∑jpj∇i,xWVj+a33∑jpj∇i,yWVj+τi,y+Fi,y,where u and v are tangential and normal velocities, respectively, and ∇i,xWand∇i,yW are kernel derivatives with respect to the particle i in x and y directions. Coefficients ααβα,β=1,2,3 are components of corrective matrix A [53]:(14)A=a11,a12,a13a21,a22,a23a31,a32,a33=∑jWVj∑jxjiWVj∑jyjiWVj∑j∇i,xWVj∑jxji∇i,xWVj∑jyji∇i,xWVj∑j∇i,yWVj∑jxji∇i,yWVj∑jyji∇i,yWVj−1.Given a well-conditioned matrixA, FPM will have higher stability and better accuracy. So, we want particle distribution to be as uniform as possible during the simulation process. Particle shifting technology (PST) can shift particles slightly through streamlines every step, thus alleviating particle clustering or sparse discretization. The PST includes the following two steps [24]:Step 1: use the following formula to modify the particle position:(15)δri=−hs22∑j1+0.2WWΔ4∇iWVj,whereδri is the displacement correction and WΔ is the value of the smoothing function in the position of the particle itself.Step 2: after completing step 1, use the following formula to correct the field variables of particles:(16)φi′=φi+δri⋅∇φi,wherei′ and i represent modified field variables and previous ones and φ stands for a general field variable.As expressed above, the FPM scheme has more advantages compared to SPH. However, FPM is usually not suitable for particles near the free surface because condition numbers here tend to be very large. At the same time, the SPH scheme naturally satisfies the dynamic boundary condition of the free-surface flows. Therefore, in this study, the SPH scheme is still applied to calculate the field variables of fluid particles near the free surface, whereas the FPM scheme is applied to the rest. Particles in the vicinity of the free surface are identified as free-surface particles and particles neighboring them, as shown in Section3.2.Here we talk about the criterion of whether a particle is located near the free surface at a certain time step [54]. First, we take the particle as the origin, make the horizontal and vertical coordinate axes, and divide its support domain into four parts (i.e., the radius of the blue circle in Figure 1 is khs), in which the number of fluid particles is Q1–Q4, respectively. Then, the coordinate axis is rotated π/4, with the number of particles in the four quadrants recorded as Q5–Q8. For the fluid particles in the vicinity of the free surface, at least one of Q1–Q8 is less than the constant K1. Otherwise, it is an interior particle. B-spline kernel function is used as the kernel function, and the corresponding constant K1 is equal to 3.Figure 1
Illustration of the method for detecting free-surface particles.
(a)(b)


## 2.4. Boundary Conditions

The particles include fluid particles and virtual particles. Virtual particles are the physical boundary in the simulation, preventing the fluid particles from passing through the boundary and providing necessary information for fluid particles in its support domain. The virtual particles involved in the simulation process can be divided into two categories according to their positions. One is virtual particles to represent the piston-type wavemaker, and the other one represents fixed solid boundary. The horizontal velocity and displacement of the virtual particles representing the wavemaker are calculated by the wave-making theory, while the position of the virtual particle forming the fixed boundary remains stable, and its velocity is determined by the velocity at the projecting point.


### 2.4.1. Solid Wall

In the following simulation, the projecting points associated with the boundary particles are created to calculate the physical information of boundary particles as shown in Figure2. The projecting point xp is determined by(17)xp=xb+εnbΔ,where xb is the position of the virtual particle relevant to the projecting point [55], nb is the normal vector of the boundary particle, and Δ is the initial particle spacing. ε is the factor to make sure that the projecting point locates at half the initial particle spacing away from the boundary line.Figure 2
Diagram of projecting points corresponding to boundary particles.The fluid variables (velocity, pressure, and density) of the projecting point can be evaluated by(18)fp=∑mj/ρjfjWpj∑mj/ρjWpj,where j represents all fluid particles in the support domain of the projecting point.Field variables of virtual particles are calculated differently from those of fluid particles. Velocity vectors are maintained by(19)ubτ=upτ,ubn=−upn,where ubτ and upτ are the tangential velocities of virtual particles and projecting points. un represents the normal velocity. Pressure and density should satisfy the weakly compressible equation as(20)pb=pp+ρ0ghb−hp,(21)ρb=ρ0+pbc2,where pb and ρb are the pressure and density of virtual particles, ppandρp are the pressure and density of projecting points, and hbandhp are the water depth of virtual particles and projecting points.In the process of numerical study, fluid particles may get close to or even pass through the boundary formed by virtual particles. To prevent this kind of error, virtual particles should exert a repulsive force on fluid particles in their support domain as(22)fij=−cui−uj∗njWijhij2njrij∗nj,ui−uj∗nj<0,0,ui−uj∗nj>0,where rij=xi−xj and hij=hsi+hsj/2, with nj the normal vector and c artificial sound speed.


### 2.4.2. Wavemaker Theory

Wave parameters are not independent of each other and should satisfy the dispersion relationship. Given water depthh, wave height H0 and wave period T are taken as variables in the following simulation process, and other wave parameters can be determined by the dispersion relationship as(23)σ2=gKtanhkh,where σ is the angular frequency, which is dictated by(24)σ=2πT.Thus, phase speedc0 and wavelength L are then available:(25)c0=gT2πtan hkh,L=c0T.The incident wave in our simulation is generated by a piston-type wavemaker located on the left side of the flume. The horizontal displacement and velocity of the piston are given, respectively [56]:(26)xpistont=s2sinσt,upistont=σs2cosσt,where s is defined by s=H0/kh. Here, H0 is the wave height, h is the water depth, and k is the wavenumber.


## 2.4.1. Solid Wall

In the following simulation, the projecting points associated with the boundary particles are created to calculate the physical information of boundary particles as shown in Figure2. The projecting point xp is determined by(17)xp=xb+εnbΔ,where xb is the position of the virtual particle relevant to the projecting point [55], nb is the normal vector of the boundary particle, and Δ is the initial particle spacing. ε is the factor to make sure that the projecting point locates at half the initial particle spacing away from the boundary line.Figure 2
Diagram of projecting points corresponding to boundary particles.The fluid variables (velocity, pressure, and density) of the projecting point can be evaluated by(18)fp=∑mj/ρjfjWpj∑mj/ρjWpj,where j represents all fluid particles in the support domain of the projecting point.Field variables of virtual particles are calculated differently from those of fluid particles. Velocity vectors are maintained by(19)ubτ=upτ,ubn=−upn,where ubτ and upτ are the tangential velocities of virtual particles and projecting points. un represents the normal velocity. Pressure and density should satisfy the weakly compressible equation as(20)pb=pp+ρ0ghb−hp,(21)ρb=ρ0+pbc2,where pb and ρb are the pressure and density of virtual particles, ppandρp are the pressure and density of projecting points, and hbandhp are the water depth of virtual particles and projecting points.In the process of numerical study, fluid particles may get close to or even pass through the boundary formed by virtual particles. To prevent this kind of error, virtual particles should exert a repulsive force on fluid particles in their support domain as(22)fij=−cui−uj∗njWijhij2njrij∗nj,ui−uj∗nj<0,0,ui−uj∗nj>0,where rij=xi−xj and hij=hsi+hsj/2, with nj the normal vector and c artificial sound speed.


## 2.4.2. Wavemaker Theory

Wave parameters are not independent of each other and should satisfy the dispersion relationship. Given water depthh, wave height H0 and wave period T are taken as variables in the following simulation process, and other wave parameters can be determined by the dispersion relationship as(23)σ2=gKtanhkh,where σ is the angular frequency, which is dictated by(24)σ=2πT.Thus, phase speedc0 and wavelength L are then available:(25)c0=gT2πtan hkh,L=c0T.The incident wave in our simulation is generated by a piston-type wavemaker located on the left side of the flume. The horizontal displacement and velocity of the piston are given, respectively [56]:(26)xpistont=s2sinσt,upistont=σs2cosσt,where s is defined by s=H0/kh. Here, H0 is the wave height, h is the water depth, and k is the wavenumber.


## 2.5. Time Integral Scheme

The Predictor-Corrector time integration algorithm is used here. At the first half step, all fluid variables are predicted using the following schemes:(27)ρin+1/2=ρin+Δt2dρindt,vin+1/2=vin+Δt2dvindt,xin+1/2=xin+Δt2vin.Then, these values at the second half step can be calculated as(28)ρ¯in+1/2=ρin+Δt2dρin+1/2dt,v¯in+1/2=vin+Δt2dvin+1/2dt,x¯in+1/2=xin+Δt2vin.Finally, all fluid variables at the next time step can be calculated as follows:(29)ρin+1=2ρ¯in+1/2−ρin,vin+1=2v¯in+1/2−vin,xin+1=2x¯in+1/2−xin.To ensure the accuracy and stability, the time step is selected by(30)dt=CFL∗Δc.Herein, the CFL number is chosen as 0.5.


## 3. Results

### 3.1. Numerical Validation: Standing Wave

Standing waves are typical free-surface flows, which have been studied by many researchers. There exists an analytic solution of total kinetic energy in its process of dissipation. In this section, we utilize coupled FPM with SPH scheme to simulate standing waves to prove its accuracy and stability. Figure3 describes the setup of our numerical model. It is a 2D rectangle flume with its top side opening and left and right side high enough. Water depth in the flume is h = Lx, where Lx is the width of the flume. At the beginning t = 0, the free surface is strictly horizontal, and the density is calculated by ρ0=ρr1+gh/c2, where c is the artificial sound speed and ρr is the reference density.Figure 3
Sketch of the validation model.Att = 0, velocity potential is given by(31)φ0x,y=−H0g2σcoshky+hcoshkhcoskx,where H0 is the wave height, σ is the wave frequency, and k is the wavenumber. Here, k is equal to 2π/L, where L is the wavelength and set to be equal to Lx in the simulation. Wave frequency is solved by the dispersion relationship. Based on equation (31), the initial velocity at every point is determined. In this case, H0 is set to be 0.1 m, while L = 1 m, c = 31.3 m/s, g = 9.8 m/s2, T = 0.801 s, and ρr = 1000.For flows without the viscosity, wave height and wave period will never change because there exists no energy dissipation. When the viscous force is considered in the simulation, the total kinetic energy decreases, resulting in less wave height until the motion stops. The kinetic energy is given by the following equation [52]:(32)Ek=−λgH0232e−βt1+cos2σt,where β = 4νk2 is the dissipation coefficient, where ν is the kinematic viscosity coefficient and ν=αc0hs/8 in 2D simulation, where hs represents smoothing length. Here, we set nx=Lx/Δ, where △ is the initial particle spacing. Timestep dt is decided by equation (30).This process is characterized by Reynolds number Re =hgh/υ. Figure 4 shows the time series of kinetic energy dissipation at Re = 5000. E0 represents kinetic energy at t = 0, which is larger than kinetic energy at any other moment throughout the process. In addition, t∗ is a dimensionless variable and is equal to t (g/h) 0.5. In Figure 4, nx = 100 means that the initial particle spacing is Lx/100 = 0.01. Figure 4 shows results predicted by coupled FPM with SPH scheme fit well with an analytic solution by (32). Obviously, as the particle spacing decreases, the results by FPM + SPH scheme gradually tend to the analytical solution, which shows good convergence.Figure 4
Time series of kinetic energy dissipation at Re = 5000.


### 3.2. Computational Model

The wave flume in the simulation is shown in Figure5. A piston wavemaker generates stable incident waves in the flume. During the operation, the function is used to monitor the wave surface changes at fixed points (G1–G13 in Figure 5). The numerical wave flume established is entirely consistent with the physical experimental model in Yao et al. [48] to study the same problem in terms of wave incident conditions (including wave parameters and model size).Figure 5
Sketch of the experimental arrangement.The size of the model is shown in Figure5. Waves are generated by the piston wavemaker. In the beginning, the wavemaker is 16.35 m away from the foot of the slope in front of the reef. The slope in front of the reef is 1 : 6 and 0.35 m high. The length of the reef flat is 20.0 m, and the slope behind the reef flat is 1 : 8. In the case of reef ridge, the reef ridge is 0.55 m long and 0.05 m high. The reason is that after passing through the breaking zone, the generated waves are very weak, the reflected waves from the far end of the flume are also very few, and the wave increase on the reef flat should be the same everywhere.As shown in G1–G13 of Figure5, the water surface height at the corresponding position is monitored during the simulation process. G1 and G2 are located 12.10 m and 12.35 m away from the pusher plate to monitor the incident wave height. G3–G5 is located at 0 m, 1 m, and 1.5 m away from the reef front slope’s toe to monitor the wave change when climbing the reef front slope. G6–G9 is located at 0 m, 0.5 m, 1 m, and 1.5 m away from the edge of reef flat, and G10–G13 is located at 1.7 m, 3.3 m, 5.0 m, and 6.7 m away from the edge of reef flat to monitor the wave height on the reef flat. There are several additional points to monitor the change of wave surface under individual conditions, which will be specially explained in the analysis. In the discussion below, if there is no particular explanation, x = 0 m refers to the horizontal position of the wavemaker at the initial time, and it is positive to the right. In the following simulation, we set Δ = 0.01 m, g = 9.8 m/s2, and ρr = 1000 at t = 0, where Δ is particle spacing and ρr is reference density unless otherwise stated. The sound speed c is set by c=10gh, where h is the water depth before the fore-reef slope.Figure6 shows the distribution of particles at the beginning with and without reef ridge. Δ is set to be 0.01 m unless otherwise stated, as is shown in the figure. The combination of different wave parameters is considered in the simulation. Water depth h selects three typical values: 0.35 m, 0.4 m, and 0.45 m. Generally, there are different submergence depths of reef flat: 0 cm, 5 cm, and 10 cm without reef ridge and 0 cm and 5 cm at the top of reef ridge with reef ridge. The study of water depth of 0.35 m in the existence of ridge is not conducted, where reef ridge separates waves and the reef flat. The incident wave period is 1.667 s, 1.25 s, and 1.0 s, and the wave height ranges from 0.046 m to 0.131 m. Due to the reflection of the slope in front of the reef, there will be a partial standing wave phenomenon in front of the slope. As a result, the actual wave height of the incident wave will be slightly different from that of the theoretical wave. In the analysis of the results, the actual incident wave height monitored in the simulation process will be used.Figure 6
Initial particle distribution for 0.45 m water depth: (a) with ridge and (b) without ridge.
(a)(b)


### 3.3. Regular Waves

As stated in Section2.3, the free-surface particles and interior particles should be separated from each other in the simulation process, or nonphysical phenomena will probably occur due to the unreasonable selection of numerical schemes. Figure 7 shows the free-surface particles distinguished by the method described in Section 2.3, in which red ones in Figure 7 are free-surface particles, gray the interior particles, and black the virtual particles. The figure is from the 200000th step of the selected typical conditions, corresponding to t = 47.62 seconds. Here, “typical conditions” mean wave depth h = 0.45 m, wave period T = 1.25 s, and wave height H0 = 0.095 m. Figure 7(a) shows the screening results of more than one full-waveform length. It can be seen that this method has high accuracy in screening free-surface particles for waves. Figure 7(b) shows a local enlargement. The 2-3 layers of particles closest to the surface will be calculated iteratively using the SPH scheme instead of the FPM scheme.Figure 7
Free-surface particles distinguished from the rest.
(a)(b)Regular waves are defined to satisfy the condition of wave heightH0 << wavelength L. Therefore, nonlinearity should be checked first before applying the theory of regular waves. According to equations (22)–(25), the wave parameters under typical conditions (h = 0.45 m, T = 1.25 s, and H0 = 0.095 m) can be obtained by using the linear wave theory (Table 1).Table 1
Wave parameters under typical conditions.
Wave parametersValueWavenumberk2.964Angular frequencyσ5.027Wave velocityc01.697 m/sWavelengthL2.122 mIt can be seen that, for typical conditions,H0/L is equal to 0.048, which does not strictly meet the linear condition. However, in the analysis of this paper, the airy wave theory is still used for the sake of convenience. Figure 8 shows the comparison results of free surface and sine function under typical working conditions (h = 0.45 m, T = 1.25 s, and H0 = 0.095 m).Figure 8
Free surface compared with a standard trigonometric functiony=0.45+0.02sin2.90x+0.03cos2.90x.Though the horizontal displacement and velocity of the piston have been proved [56], the expected wave parameters set on the wavemaker are slightly different from the incident wave parameters generated actually in the simulation (i.e., the incident wave parameters obtained by monitoring the wave conditions at G1 and G2 positions). One of the essential reasons is that the waves will be reflected between the wavemaker and the fore-reef slope, possibly forming a particular standing wave phenomenon under some conditions. Therefore, the actual incident wave parameters should not be deduced from the preset wave parameters on the wavemaker before simulation. Instead, the recorded waves measured at G1 and G2 should be used for separating the incident wave from the reflected wave to obtain the actual incident wave, which is used to analyze the subsequent wave changes.Figure9 presents the time series of wave records at G1 (x = 12.10 m) and G2 (x = 12.35 m) in typical conditions. It can be seen that, in the wave process simulated, the water surface duration curves of G1 and G2 are different only in their phase, while the wave height is the same, so it is assumed that there is almost no reflection, and the influence of wave reflection on the incident wave is not considered in the following analysis.Figure 9
Time series of wave records at G1 and G2.


### 3.4. Wave Breaking

Generally, there are three main forms of wave breaking: spilling, plunging, and surging. In the physical model test, the type of breaking and the beginning time of breaking can be recognized easily: the spilling wave breaking begins when the white breaker appears at the wave crest, and the corresponding wave height is the breaking wave height. For the plunging wave breaking, it is considered that the breaking starts when the wave crest is near vertical. For the surging wave breaking, the appearance of spray at the lower part of the wavefront peak is considered as the beginning of breaking. However, in our simulation, the breaking types can only be classified by the change of wave shape: under most tested conditions, the breaking of waves is in the form of plunging breaking, and another form of breaking is observed when the water depth in front of some reefs is deeper (water depthh = 0.45 m) with a small wave height. It can be judged that it should belong to a spilling breaker.Figure10 shows the wave propagating and breaking process near the fore-reef slope in front of the reef (h = 0.40 m, T = 1.667 s, and H0 = 0.107 m) and the pressure distribution at the corresponding time. According to equation (3), the pressure has a linear relationship with density. Thus, the fluid density distribution can also be seen. The wave breaking process has several apparent stages. In Figure 10(a), due to the rapid decrease of water depth and the increase of nonlinear influence, the two sides of the wave crest become extremely asymmetric, but there is no breaker at this time. In Figure 10(b), the peak reverses to form a tongue shape with aeration, and the pressure and density are extremely uneven. In Figure 10(c), the tongue-shaped wave crest submerges into the water, forming strong turbulence and moving forward. In Figure 10(d), the waves reach the shoreward side of the surf zone, gradually generating new waves. However, the waves at this time are still quite nonlinear until outside the surf zone.Figure 10
The transforming and breaking of plunging breakers (h = 0.40 m, t = 1.667 s, and H0 = 0.107 m).Figure11 gives time series of wave records at different positions under the same wave condition (h = 0.45 m, T = 1.25 s, and H0 = 0.095 m). It can be seen that, on the seaward side of the slope in front of the reef, its shape is close to regular waves when the wave has not climbed the slope. The wave gradually becomes nonlinear after climbing the slope for a certain distance. In the case of reef ridge, the wave height here is smaller than that without reef ridge. The possible reason is that the vertical surface of the ridge increases the wave reflection and thus interferes with the incident wave at this point. In Figure 11(c), G7 is in the area with the strongest turbulence during the wave propagation in the whole reef terrain: for the terrain without reef ridge, the tongue of the breaking wave is formed nearby, and the air is mixed into the water. The wave has just undergone a sudden change of water depth for the terrain with reef ridge, and the slope quickly rushes up to the reef ridge. In this region, the detected waveform is not even the waveform of solitary wave or Stokes wave, but a kind of sawtooth. In Figure 11(d), G10 is already at the end of the surf zone. After the strong turbulence in the surf zone, the waves gradually took on a conical wave shape. After this point, the wave will be adjusted to a regular waveform. But at the same time, the wave height is much smaller than the incident wave, and the wavelength is longer because of the crushing dissipation of a large amount of energy.Figure 11
Time series of wave records at different locations. (a) G1. (b) G4. (c) G7. (d) G10.
(a)(b)(c)(d)


### 3.5. Wave-Induced Setup

In the simulation process, 13 gauges (G1–G13) in the flume can record the wave surface at a fixed time interval. Figure12 shows the average water surface height of each point on the reef when h = 0.45 m, T = 1.25 s, and H0 = 0.095 m. In the figure, x = 0 is the position of the slope toe in front of the reef, and the shoreward is positive. In addition, the experimental results of Yao et al. [17] on the same problem are given for comparison. No matter whether there is a reef ridge or not, the trend of setup and setdown is consistent. Setdown occurs at the seaward side of the surf zone. After that, waves arrive at the surf zone, where the mean water level rises rapidly until stability and remains unchanged after 5 m away from the slope toe. When the reef ridge exists, there will be a larger setup on the reef flat. However, numerical simulation results are quite different from the experimental results within 0–2 m from the slope toe, and the law reflected is also inconsistent with the existing theory. The reasons will be discussed in the following.Figure 12
Mean water level (MWL) offshore and across the reef profile under different wave conditions.Figure13 shows the comparison of setup amplitude of reef flat without reef ridge under different wave conditions, including various water depths in front of the reef (0.35 m, 0.40 m, and 0.45 m), incident wave periods (1.0 s, 1.25 s, and 1.667 s), and incident wave heights (0.054 m–0.131 m). It can be seen that the waves have a very significant effect on the setup of reef flat: there is a positive correlation between the wave height of the incident wave and the water depth in front of the reef. Particularly, the relationship between the setup amplitude and the incident wave height is even close to linear. On the other hand, the increase or decrease of the water depth in front of the reef will affect the effects of the other two factors: when the water depth in front of the reef is small, with the increase of wave period and wave height, the setup of the reef flat increases significantly. With the increase of water depth, the influence of wave period and wave height on setup becomes less significant. In addition, experimental results when black markers give h = 0.4 m by Yao et al. [17] in this figure to compare, which correspond with numerical results in the same condition, reveal that the SPH-FPM method is capable of dealing with wave dynamics.Figure 13
Maximum wave setup on reef flat as a function of deep-water wave height for different wave periods, water depths, and incident wave heights.Figure14 shows the comparison of the effects of reef ridge on the setup of reef flat under different wave conditions with a 0.45 m water depth in front of the reef. It can be seen that the existence of reef ridge can significantly increase the setup of the reef flat.Figure 14
Variation of wave setup with the existence of reef ridge.


## 3.1. Numerical Validation: Standing Wave

Standing waves are typical free-surface flows, which have been studied by many researchers. There exists an analytic solution of total kinetic energy in its process of dissipation. In this section, we utilize coupled FPM with SPH scheme to simulate standing waves to prove its accuracy and stability. Figure3 describes the setup of our numerical model. It is a 2D rectangle flume with its top side opening and left and right side high enough. Water depth in the flume is h = Lx, where Lx is the width of the flume. At the beginning t = 0, the free surface is strictly horizontal, and the density is calculated by ρ0=ρr1+gh/c2, where c is the artificial sound speed and ρr is the reference density.Figure 3
Sketch of the validation model.Att = 0, velocity potential is given by(31)φ0x,y=−H0g2σcoshky+hcoshkhcoskx,where H0 is the wave height, σ is the wave frequency, and k is the wavenumber. Here, k is equal to 2π/L, where L is the wavelength and set to be equal to Lx in the simulation. Wave frequency is solved by the dispersion relationship. Based on equation (31), the initial velocity at every point is determined. In this case, H0 is set to be 0.1 m, while L = 1 m, c = 31.3 m/s, g = 9.8 m/s2, T = 0.801 s, and ρr = 1000.For flows without the viscosity, wave height and wave period will never change because there exists no energy dissipation. When the viscous force is considered in the simulation, the total kinetic energy decreases, resulting in less wave height until the motion stops. The kinetic energy is given by the following equation [52]:(32)Ek=−λgH0232e−βt1+cos2σt,where β = 4νk2 is the dissipation coefficient, where ν is the kinematic viscosity coefficient and ν=αc0hs/8 in 2D simulation, where hs represents smoothing length. Here, we set nx=Lx/Δ, where △ is the initial particle spacing. Timestep dt is decided by equation (30).This process is characterized by Reynolds number Re =hgh/υ. Figure 4 shows the time series of kinetic energy dissipation at Re = 5000. E0 represents kinetic energy at t = 0, which is larger than kinetic energy at any other moment throughout the process. In addition, t∗ is a dimensionless variable and is equal to t (g/h) 0.5. In Figure 4, nx = 100 means that the initial particle spacing is Lx/100 = 0.01. Figure 4 shows results predicted by coupled FPM with SPH scheme fit well with an analytic solution by (32). Obviously, as the particle spacing decreases, the results by FPM + SPH scheme gradually tend to the analytical solution, which shows good convergence.Figure 4
Time series of kinetic energy dissipation at Re = 5000.


## 3.2. Computational Model

The wave flume in the simulation is shown in Figure5. A piston wavemaker generates stable incident waves in the flume. During the operation, the function is used to monitor the wave surface changes at fixed points (G1–G13 in Figure 5). The numerical wave flume established is entirely consistent with the physical experimental model in Yao et al. [48] to study the same problem in terms of wave incident conditions (including wave parameters and model size).Figure 5
Sketch of the experimental arrangement.The size of the model is shown in Figure5. Waves are generated by the piston wavemaker. In the beginning, the wavemaker is 16.35 m away from the foot of the slope in front of the reef. The slope in front of the reef is 1 : 6 and 0.35 m high. The length of the reef flat is 20.0 m, and the slope behind the reef flat is 1 : 8. In the case of reef ridge, the reef ridge is 0.55 m long and 0.05 m high. The reason is that after passing through the breaking zone, the generated waves are very weak, the reflected waves from the far end of the flume are also very few, and the wave increase on the reef flat should be the same everywhere.As shown in G1–G13 of Figure5, the water surface height at the corresponding position is monitored during the simulation process. G1 and G2 are located 12.10 m and 12.35 m away from the pusher plate to monitor the incident wave height. G3–G5 is located at 0 m, 1 m, and 1.5 m away from the reef front slope’s toe to monitor the wave change when climbing the reef front slope. G6–G9 is located at 0 m, 0.5 m, 1 m, and 1.5 m away from the edge of reef flat, and G10–G13 is located at 1.7 m, 3.3 m, 5.0 m, and 6.7 m away from the edge of reef flat to monitor the wave height on the reef flat. There are several additional points to monitor the change of wave surface under individual conditions, which will be specially explained in the analysis. In the discussion below, if there is no particular explanation, x = 0 m refers to the horizontal position of the wavemaker at the initial time, and it is positive to the right. In the following simulation, we set Δ = 0.01 m, g = 9.8 m/s2, and ρr = 1000 at t = 0, where Δ is particle spacing and ρr is reference density unless otherwise stated. The sound speed c is set by c=10gh, where h is the water depth before the fore-reef slope.Figure6 shows the distribution of particles at the beginning with and without reef ridge. Δ is set to be 0.01 m unless otherwise stated, as is shown in the figure. The combination of different wave parameters is considered in the simulation. Water depth h selects three typical values: 0.35 m, 0.4 m, and 0.45 m. Generally, there are different submergence depths of reef flat: 0 cm, 5 cm, and 10 cm without reef ridge and 0 cm and 5 cm at the top of reef ridge with reef ridge. The study of water depth of 0.35 m in the existence of ridge is not conducted, where reef ridge separates waves and the reef flat. The incident wave period is 1.667 s, 1.25 s, and 1.0 s, and the wave height ranges from 0.046 m to 0.131 m. Due to the reflection of the slope in front of the reef, there will be a partial standing wave phenomenon in front of the slope. As a result, the actual wave height of the incident wave will be slightly different from that of the theoretical wave. In the analysis of the results, the actual incident wave height monitored in the simulation process will be used.Figure 6
Initial particle distribution for 0.45 m water depth: (a) with ridge and (b) without ridge.
(a)(b)


## 3.3. Regular Waves

As stated in Section2.3, the free-surface particles and interior particles should be separated from each other in the simulation process, or nonphysical phenomena will probably occur due to the unreasonable selection of numerical schemes. Figure 7 shows the free-surface particles distinguished by the method described in Section 2.3, in which red ones in Figure 7 are free-surface particles, gray the interior particles, and black the virtual particles. The figure is from the 200000th step of the selected typical conditions, corresponding to t = 47.62 seconds. Here, “typical conditions” mean wave depth h = 0.45 m, wave period T = 1.25 s, and wave height H0 = 0.095 m. Figure 7(a) shows the screening results of more than one full-waveform length. It can be seen that this method has high accuracy in screening free-surface particles for waves. Figure 7(b) shows a local enlargement. The 2-3 layers of particles closest to the surface will be calculated iteratively using the SPH scheme instead of the FPM scheme.Figure 7
Free-surface particles distinguished from the rest.
(a)(b)Regular waves are defined to satisfy the condition of wave heightH0 << wavelength L. Therefore, nonlinearity should be checked first before applying the theory of regular waves. According to equations (22)–(25), the wave parameters under typical conditions (h = 0.45 m, T = 1.25 s, and H0 = 0.095 m) can be obtained by using the linear wave theory (Table 1).Table 1
Wave parameters under typical conditions.
Wave parametersValueWavenumberk2.964Angular frequencyσ5.027Wave velocityc01.697 m/sWavelengthL2.122 mIt can be seen that, for typical conditions,H0/L is equal to 0.048, which does not strictly meet the linear condition. However, in the analysis of this paper, the airy wave theory is still used for the sake of convenience. Figure 8 shows the comparison results of free surface and sine function under typical working conditions (h = 0.45 m, T = 1.25 s, and H0 = 0.095 m).Figure 8
Free surface compared with a standard trigonometric functiony=0.45+0.02sin2.90x+0.03cos2.90x.Though the horizontal displacement and velocity of the piston have been proved [56], the expected wave parameters set on the wavemaker are slightly different from the incident wave parameters generated actually in the simulation (i.e., the incident wave parameters obtained by monitoring the wave conditions at G1 and G2 positions). One of the essential reasons is that the waves will be reflected between the wavemaker and the fore-reef slope, possibly forming a particular standing wave phenomenon under some conditions. Therefore, the actual incident wave parameters should not be deduced from the preset wave parameters on the wavemaker before simulation. Instead, the recorded waves measured at G1 and G2 should be used for separating the incident wave from the reflected wave to obtain the actual incident wave, which is used to analyze the subsequent wave changes.Figure9 presents the time series of wave records at G1 (x = 12.10 m) and G2 (x = 12.35 m) in typical conditions. It can be seen that, in the wave process simulated, the water surface duration curves of G1 and G2 are different only in their phase, while the wave height is the same, so it is assumed that there is almost no reflection, and the influence of wave reflection on the incident wave is not considered in the following analysis.Figure 9
Time series of wave records at G1 and G2.


## 3.4. Wave Breaking

Generally, there are three main forms of wave breaking: spilling, plunging, and surging. In the physical model test, the type of breaking and the beginning time of breaking can be recognized easily: the spilling wave breaking begins when the white breaker appears at the wave crest, and the corresponding wave height is the breaking wave height. For the plunging wave breaking, it is considered that the breaking starts when the wave crest is near vertical. For the surging wave breaking, the appearance of spray at the lower part of the wavefront peak is considered as the beginning of breaking. However, in our simulation, the breaking types can only be classified by the change of wave shape: under most tested conditions, the breaking of waves is in the form of plunging breaking, and another form of breaking is observed when the water depth in front of some reefs is deeper (water depthh = 0.45 m) with a small wave height. It can be judged that it should belong to a spilling breaker.Figure10 shows the wave propagating and breaking process near the fore-reef slope in front of the reef (h = 0.40 m, T = 1.667 s, and H0 = 0.107 m) and the pressure distribution at the corresponding time. According to equation (3), the pressure has a linear relationship with density. Thus, the fluid density distribution can also be seen. The wave breaking process has several apparent stages. In Figure 10(a), due to the rapid decrease of water depth and the increase of nonlinear influence, the two sides of the wave crest become extremely asymmetric, but there is no breaker at this time. In Figure 10(b), the peak reverses to form a tongue shape with aeration, and the pressure and density are extremely uneven. In Figure 10(c), the tongue-shaped wave crest submerges into the water, forming strong turbulence and moving forward. In Figure 10(d), the waves reach the shoreward side of the surf zone, gradually generating new waves. However, the waves at this time are still quite nonlinear until outside the surf zone.Figure 10
The transforming and breaking of plunging breakers (h = 0.40 m, t = 1.667 s, and H0 = 0.107 m).Figure11 gives time series of wave records at different positions under the same wave condition (h = 0.45 m, T = 1.25 s, and H0 = 0.095 m). It can be seen that, on the seaward side of the slope in front of the reef, its shape is close to regular waves when the wave has not climbed the slope. The wave gradually becomes nonlinear after climbing the slope for a certain distance. In the case of reef ridge, the wave height here is smaller than that without reef ridge. The possible reason is that the vertical surface of the ridge increases the wave reflection and thus interferes with the incident wave at this point. In Figure 11(c), G7 is in the area with the strongest turbulence during the wave propagation in the whole reef terrain: for the terrain without reef ridge, the tongue of the breaking wave is formed nearby, and the air is mixed into the water. The wave has just undergone a sudden change of water depth for the terrain with reef ridge, and the slope quickly rushes up to the reef ridge. In this region, the detected waveform is not even the waveform of solitary wave or Stokes wave, but a kind of sawtooth. In Figure 11(d), G10 is already at the end of the surf zone. After the strong turbulence in the surf zone, the waves gradually took on a conical wave shape. After this point, the wave will be adjusted to a regular waveform. But at the same time, the wave height is much smaller than the incident wave, and the wavelength is longer because of the crushing dissipation of a large amount of energy.Figure 11
Time series of wave records at different locations. (a) G1. (b) G4. (c) G7. (d) G10.
(a)(b)(c)(d)


## 3.5. Wave-Induced Setup

In the simulation process, 13 gauges (G1–G13) in the flume can record the wave surface at a fixed time interval. Figure12 shows the average water surface height of each point on the reef when h = 0.45 m, T = 1.25 s, and H0 = 0.095 m. In the figure, x = 0 is the position of the slope toe in front of the reef, and the shoreward is positive. In addition, the experimental results of Yao et al. [17] on the same problem are given for comparison. No matter whether there is a reef ridge or not, the trend of setup and setdown is consistent. Setdown occurs at the seaward side of the surf zone. After that, waves arrive at the surf zone, where the mean water level rises rapidly until stability and remains unchanged after 5 m away from the slope toe. When the reef ridge exists, there will be a larger setup on the reef flat. However, numerical simulation results are quite different from the experimental results within 0–2 m from the slope toe, and the law reflected is also inconsistent with the existing theory. The reasons will be discussed in the following.Figure 12
Mean water level (MWL) offshore and across the reef profile under different wave conditions.Figure13 shows the comparison of setup amplitude of reef flat without reef ridge under different wave conditions, including various water depths in front of the reef (0.35 m, 0.40 m, and 0.45 m), incident wave periods (1.0 s, 1.25 s, and 1.667 s), and incident wave heights (0.054 m–0.131 m). It can be seen that the waves have a very significant effect on the setup of reef flat: there is a positive correlation between the wave height of the incident wave and the water depth in front of the reef. Particularly, the relationship between the setup amplitude and the incident wave height is even close to linear. On the other hand, the increase or decrease of the water depth in front of the reef will affect the effects of the other two factors: when the water depth in front of the reef is small, with the increase of wave period and wave height, the setup of the reef flat increases significantly. With the increase of water depth, the influence of wave period and wave height on setup becomes less significant. In addition, experimental results when black markers give h = 0.4 m by Yao et al. [17] in this figure to compare, which correspond with numerical results in the same condition, reveal that the SPH-FPM method is capable of dealing with wave dynamics.Figure 13
Maximum wave setup on reef flat as a function of deep-water wave height for different wave periods, water depths, and incident wave heights.Figure14 shows the comparison of the effects of reef ridge on the setup of reef flat under different wave conditions with a 0.45 m water depth in front of the reef. It can be seen that the existence of reef ridge can significantly increase the setup of the reef flat.Figure 14
Variation of wave setup with the existence of reef ridge.


## 4. Discussion

This section discusses the effects of different wave parameters on the experimental results. In the simulation, due to the conservation of mass in the closed flume, the average water level of the front reef dropped (see Figure12). However, this decline is generally less than 1 cm. The typical test conditions with a wave period of 1.25 s and a water depth of 45 cm will only reduce the wavelength by about 1% according to the dispersion equation, so the impact of water level reduction on the characteristics of incident wave can be ignored. It can also be seen that when comparing the increase and decrease of water level caused by waves, it is reasonable to define them as the deviation between the mean water level and the still water level without considering the drop of water level in front of reefs.As shown in Figure12, the simulation results are partly different from those in the experiment by Yao et al. [48], especially in the part of wave-induced setdown. It is probably because of the following:(1)
In most conditions, the generated wave is close to linear waves, but it does not strictly meet the linear condition that the wave heightH0 is far less than the wavelength L. The actual H0/L = 1/10∼1/30, but the linear wave theory is still used instead of Stokes wave theory.(2)
As shown in Figure9, the reflection of incident waves on the slope in front of the reef has little interference in itself. However, this interference should be considered, especially in the case of reef ridge. When the ridge exists, the measured wave reflection coefficient Kr is between 6% and 55% [48]. That is to say, the two-point method should be used to separate the actual incident wave from the monitored wave surface curve in the conditions with a large reflection coefficient. Otherwise, there will be large errors. This may also explain that the calculated maximum wave-induced setdown without reef ridge is closer to the experimental result, while it deviates greatly from the experimental value when the ridge exists.(3)
Due to limited calculation time and too many calculated conditions, the simulation is only limited to 0–45 s since the beginning. Taking the typical condition as an example, the wave velocityc0 is 1.7 m/s. Thus, the wave propagates about 2.4 rounds in front of the reef (16.35 m one way) during this period. However, Yao et al. [48] found that the waves reach a stable state after about 3 min since the beginning of wave generation. Although the wave has been reflected several times in front of the reef and the shape of the wave has gradually been stabilized, other processes, especially the water mass transportation from the front of the reef to the reef flat, might be far from completeness. In Figure 12, an obvious water setdown can be observed at x = 2. But at x = 0, far away from the reef flat, the numerical result is not as good as the experimental result, probably because the mass transportation process has not been completed.


## 5. Conclusions

Thanks to their particular topography, reefs can greatly reduce the impact of hazards like storm surge to protect the coastal areas. In this study, the improved SPH method is used to simulate the process of 2D regular waves transforming and breaking on reefs. In the simulation, the parameters including water depth in front of the fore-reef (submergence over the reef flat), wave period, wave height, and the presence of reef ridge are compared with the experimental results and discussed. The simulation results show the following:(1)
The wave propagation over the reef experiences the process of shoaling, breaking, and reforming. Wave breaking is mainly in the form of plunging breaking with obvious stages: when climbing the slope, the two sides of the wave ridge become extremely asymmetric until the wave crest turns back, forming a tongue shape; after continuous propagation, the tongue-shaped wave crest submerges into the water, resulting in strong turbulence; finally, the waves will reach the shoreward side of the surf zone, where they gradually begin to regenerate new waves.(2)
Wave breaking is accompanied by setup after the reef ridge and setdown before it. No matter whether there is a reef ridge or not, the trend of wave-induced setup is similar: the mean water level decreases slightly when the waves begin to climb the fore-reef slope. After that, the mean water level rises rapidly until it reaches a stable range and remains unchanged after 5 m away from the slope toe.(3)
With the increase of incident wave period and wave height, as well as the decrease of water depth in front of the reef, the setup on the reef flat will increase correspondingly. The existence of reef ridge will make the wave breaking more intense, dissipate more energy in the breaking process, and thus increase the setup.



---

*Source: 9669925-2022-05-05.xml*