Characteristics of Breaking Wave Forces on Piles over a Permeable Seabed

: Most offshore wind turbines are installed in shallow water and exposed to breaking waves. Previous numerical studies focusing on breaking wave forces generally ignored the seabed permeability. In this paper, a numerical model based on Volume-Averaged Reynolds Averaged Navier–Stokes equations (VARANS) is employed to reveal the process of a solitary wave interacting with a rigid pile over a permeable slope. Through applying the Forchheimer saturated drag equation, effects of seabed permeability on ﬂuid motions are simulated. The reliability of the present model is veriﬁed by comparisons between experimentally obtained data and the numerical results. Further, 190 cases are simulated and the effects of different parameters on breaking wave forces on the pile are studied systematically. Results indicate that over a permeable seabed, the maximum breaking wave forces can occur not only when waves break just before the pile, but also when a “secondary wave wall” slams against the pile, after wave breaking. With the initial wave height increasing, breaking wave forces will increase, but the growth can decrease as the slope angle and permeability increase. For inclined piles around the wave breaking point, the maximum breaking wave force usually occurs with an inclination angle of α = − 22.5 ◦ or 0 ◦ . of a solitary wave permeable


Introduction
Most offshore wind turbines (OWTs) are installed in shallow-water zones at a depth of less than 30 m [1] and subject to the impact of breaking waves. The breaking waves slamming against marine structures can impose an additional wave force on them, at a maximum of 8 to 9 times than the force caused by non-breaking waves [2], threatening the structural safety of OWT foundations.
To date, the pile is the most popular foundation form of OWTs in Europe. WindEurope's statistics show that 80% of European OWTs have adopted pile foundations [3]. While in China, the high-rise pile cap foundation, which can be regarded as a combination of several inclined piles and a concrete pile cap, is utilized more and more frequently due to its great stability and economic efficiency in soft clay ( Figure 1). Hence, except for the interaction between seabed soil and pile [4][5][6], breaking wave forces acting on vertical and inclined piles have always been a major concern in the field of science and engineering [7][8][9][10][11][12][13][14][15][16][17][18][19].
With the increase in computational capabilities, the high-precision computational fluid dynamics (CFD) simulation provides an important means to investigate the evolution of wave shoaling and breaking. Mo et al. [7] reported both numerical and experimental investigations into the runup process of solitary waves and the time histories of wave forces acting on a cylinder. Based on Reynolds Averaged Navier-Stokes equations (RANS) with a k-ε model for turbulence closure, Xiao et al. [8] built a three-dimensional (3-D) numerical flume and investigated the variations in breaking wave forces with relative positions of the pile to the shoreline. Choi et al. [9] came up with a method for obtaining net breaking wave impact forces on vertical and inclined piles from the measured response time histories, based on the modified Navier-Stokes (N-S) equations. Kamath et al. [10] and Chella et al. [11] took the effect of wave shapes into account, and studied the characteristics of breaking wave forces when the relative position of piles to the breaking point changed, utilizing the open source software REEF3D. Further, Chella et al. [12] revealed the influence of the incident wave characteristics, wave impact conditions and breaker types on the maximum wave impact pressure and kinematics. Besides, Johlas et al. [13] evaluated and improved the accuracy of common breaking wave criteria through comparison to CFD simulations of breaking waves. Lomonaco et al. [14] provided a detailed data set for developing empirical predictive models of wave loads on offshore wind piles. J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 2 of 18 wave forces acting on a cylinder. Based on Reynolds Averaged Navier-Stokes equations (RANS) with a k-ε model for turbulence closure, Xiao et al. [8] built a three-dimensional (3-D) numerical flume and investigated the variations in breaking wave forces with relative positions of the pile to the shoreline. Choi et al. [9] came up with a method for obtaining net breaking wave impact forces on vertical and inclined piles from the measured response time histories, based on the modified Navier-Stokes (N-S) equations. Kamath et al. [10] and Chella et al. [11] took the effect of wave shapes into account, and studied the characteristics of breaking wave forces when the relative position of piles to the breaking point changed, utilizing the open source software REEF3D. Further, Chella et al. [12] revealed the influence of the incident wave characteristics, wave impact conditions and breaker types on the maximum wave impact pressure and kinematics. Besides, Johlas et al. [13] evaluated and improved the accuracy of common breaking wave criteria through comparison to CFD simulations of breaking waves. Lomonaco et al. [14] provided a detailed data set for developing empirical predictive models of wave loads on offshore wind piles. All studies above assume the seabed to be impermeable. However, the real seabed is made of porous media, and its permeability can affect not only the process of wave breaking [20][21][22] but also the breaking wave force. Huang et al. [23] simulated the propagation of a solitary wave over a permeable seabed using a two-dimensional (2-D) N-S model. Hammeken [24] presented experimental and numerical studies to investigate the swash hydrodynamics and their interaction with a permeable seabed, and found that the introduction of seabed permeability could decrease the distance between the breaking points and the shoreline for plunging waves.
However, most previous studies related to permeable seabed focused on porous breakwaters or breakwaters on rubble mound [25][26][27], little attention has been paid to the piles. Chen et al. [28] investigated the effect of seabed permeability on wave forces over a flat seabed, whereas the wave breaking was not taken into consideration.
Since the permeability of a porous slope has not been considered in previous literature related to wave breaking, the knowledge on the effects of some key factors, such as seabed permeability, incident wave heights and slope angles, on breaking wave forces acting on a pile over a permeable slope is still missing. Hence, the investigation on the effects of seabed permeability on breaking wave forces makes significant sense.
In this paper, a 3-D numerical model based on Volume-Averaged Reynolds Averaged Navier-Stokes (VARANS) equations with a k-ε model for turbulence closure is utilized. Through applying the Forchheimer saturated drag equation, the effects of seabed All studies above assume the seabed to be impermeable. However, the real seabed is made of porous media, and its permeability can affect not only the process of wave breaking [20][21][22] but also the breaking wave force. Huang et al. [23] simulated the propagation of a solitary wave over a permeable seabed using a two-dimensional (2-D) N-S model. Hammeken [24] presented experimental and numerical studies to investigate the swash hydrodynamics and their interaction with a permeable seabed, and found that the introduction of seabed permeability could decrease the distance between the breaking points and the shoreline for plunging waves.
However, most previous studies related to permeable seabed focused on porous breakwaters or breakwaters on rubble mound [25][26][27], little attention has been paid to the piles. Chen et al. [28] investigated the effect of seabed permeability on wave forces over a flat seabed, whereas the wave breaking was not taken into consideration.
Since the permeability of a porous slope has not been considered in previous literature related to wave breaking, the knowledge on the effects of some key factors, such as seabed permeability, incident wave heights and slope angles, on breaking wave forces acting on a pile over a permeable slope is still missing. Hence, the investigation on the effects of seabed permeability on breaking wave forces makes significant sense.
In this paper, a 3-D numerical model based on Volume-Averaged Reynolds Averaged Navier-Stokes (VARANS) equations with a k-ε model for turbulence closure is utilized. Through applying the Forchheimer saturated drag equation, the effects of seabed permeability on fluid motions are simulated and assessed. The numerical model is validated by comparing the calculated free surface elevations, velocity distributions and wave forces with experimental data. Following validation, 190 cases are simulated to reveal the process of a solitary wave propagating on a permeable slope and interacting with a rigid pile.
The effects of seabed permeability on breaking wave forces exerted on a pile have been investigated systematically. The influence of different input parameters on breaking wave forces acting on a pile founded over a permeable seabed is also studied and discussed.

Numerical Model Implementation
The fluid flow problem is solved with the VARANS equations. Turbulence modeling is handled using the two-equation k-ε model, which is adopted widely because of its accuracy and efficiency [8,27,29,30]. The Volume-of-Fluid (VOF) method is used to capture the free surface. With introducing the Forchheimer equation, the hindering effect on fluid caused by the permeable seabed is simulated.

Momentum equation
∂u ∂t ∂v ∂t ∂w ∂t where (x, y, z) are the axes in the Cartesian coordinates; (u, v, w) are the velocity components in the x, y, and z directions, respectively; p is water pressure; ρ is water density; V F is the fractional volume open to flow; (A x , A y , A z ) are the fractional areas open to flow in the x, y, and z directions, respectively; (G x , G y , G z ) are the components of body accelerations; (f x , f y , f z ) are the components of viscous accelerations; and (b x , b y , b z ) are the components of the drag term caused by the porous media [30,31]. The standard two-equation k-ε model [32] is used for turbulence closure: ∂k ∂t ∂ε ∂t where k is the turbulent kinetic energy; ε is its dissipation; M is the buoyancy production term, whose value is 0.0; C 1ε , C 2ε , and C 3ε are dimensionless parameters, and have defaults of 1.44, 1.92 and 0.2 [32]. The turbulent kinetic energy production N is expressed as: where C sp is a turbulence parameter, whose default value is 1.0; µ is the molecular dynamic viscosity. D k and D ε are the diffusion terms caused by the turbulent kinetic energy and its dissipation, respectively: where v k is the diffusion coefficient of k, computed based on the local value of the turbulent viscosity.

Free Surface Capture Method
The VOF method is introduced to capture the free surface. With the calculated fractional volume of fluid in every control volume, the position of the free surface is obtained: if V F = 0, the control volume is empty; if V F = 1, the control volume is full of fluid; and if 0 < V F < 1, the control volume must contain the free surface. The continuity equation of the fractional volume is as follows [33]: where u i is the velocity vector of mean flow; x i is the coordinate vector.

Permeable Seabed Model
In order to represent the effect of seabed permeability on fluid motions, the drag term b i , which is proportional to the fluid velocity, is added to the right side of the momentum Equations (2)-(4) [34]: where F d is the porous media drag coefficient; U m is the pore velocity: where U is the seepage flow velocity; φ is the porosity of the porous media. Broken waves slam on the beach slope directly on the permeable seabed, leading to a relatively high flow velocity inside the porous media. The flow inside the porous media can no longer be described by Darcy's law and an alternative porous media theory should be introduced. Thus, the drag coefficient F d is calculated using the Forchheimer saturated drag model [34]. Previous studies show that this model provides sufficient accuracy required to simulate the drag force caused by the porous media [35][36][37][38]: where R ep is the pore Reynolds number; d pore is the average pore throat diameter within the porous media; A and B are the linear and the non-linear drag coefficients, related to the experimentally determined coefficients a and b by the following equations [39]: When test data is unavailable, A and B can be estimated by [39]: where α is a constant, typically of order 180; β is a roughness factor typically ranging between 1.8 and 4.0, and for this study its value is 3.0, representing moderate surface roughness of the seabed; d 50 is the median particle diameter of the porous media. To derive the permeability K from the geometrical parameters of the porous media, the Kozeny-Carman equation proposed by Carman [40] is used: Mo et al. [7] performed a flume experiment to investigate the evolution of solitary waves breaking and interacting with a vertical cylinder resting on an impermeable slope. Particle Image Velocimetry (PIV) measurements were conducted to record wave shapes and velocity profiles. Based on the present model, a numerical case which replicated the experimental configuration in Mo et al. [7] was simulated. The grid size was 0.005 m, and the initial time step was 0.005 s. Figure 2 shows a comparison of the results calculated by the model presented here versus the experimental data in Mo et al. [7] at two different instants. Note that the time t is normalized by h 0 /g , and h 0 is the initial water depth.
where α is a constant, typically of order 180; β is a roughness factor typically ranging between 1.8 and 4.0, and for this study its value is 3.0, representing moderate surface roughness of the seabed; d50 is the median particle diameter of the porous media. To derive the permeability K from the geometrical parameters of the porous media, the Kozeny-Carman equation proposed by Carman [40] is used:

Comparison with Experimental Data on Free Surface Elevations and Velocities of Waves Propagating on an Impermeable Slope
Mo et al. [7] performed a flume experiment to investigate the evolution of solitary waves breaking and interacting with a vertical cylinder resting on an impermeable slope. Particle Image Velocimetry (PIV) measurements were conducted to record wave shapes and velocity profiles. Based on the present model, a numerical case which replicated the experimental configuration in Mo et al. [7] was simulated. The grid size was 0.005 m, and the initial time step was 0.005 s. Figure 2 shows a comparison of the results calculated by the model presented here versus the experimental data in Mo et al. [7] at two different instants. Note that the time t is normalized by 0 / h g , and h0 is the initial water depth.  [7]. Reproduced from Mo et al. [7]; with permission from Elsevier, 2021.
As Figure 2 indicates, while waves are propagating, breaking, and slamming against the pile, the present model can simulate the free surface elevations precisely, and gain relative accuracy to capture the vertical distribution of the horizontal velocities at specific positions. The calculated development of the velocities near the flume bottom is also reasonably satisfactory. It is only near the free surface where the calculated velocities are larger than the experimental data, reflecting the evident characteristic of wave breaking.

Comparison with Numerical Results on Breaking Wave Forces Acting on a Vertical Cylinder
Chella et al. [11] simulated breaking wave forces at different cylinder positions relative to the breaking point on an impermeable slope, and analyzed the effects of different wave shapes on the wave forces. The normalized relative distance Lc, used to distinguish the different relative positions of cylinders, is defined as:  [7]. Reproduced from Mo et al. [7]; with permission from Elsevier, 2021.
As Figure 2 indicates, while waves are propagating, breaking, and slamming against the pile, the present model can simulate the free surface elevations precisely, and gain relative accuracy to capture the vertical distribution of the horizontal velocities at specific positions. The calculated development of the velocities near the flume bottom is also reasonably satisfactory. It is only near the free surface where the calculated velocities are larger than the experimental data, reflecting the evident characteristic of wave breaking.

Comparison with Numerical Results on Breaking Wave Forces Acting on a Vertical Cylinder
Chella et al. [11] simulated breaking wave forces at different cylinder positions relative to the breaking point on an impermeable slope, and analyzed the effects of different wave shapes on the wave forces. The normalized relative distance L c , used to distinguish the different relative positions of cylinders, is defined as: where H b is the breaker height; x c is the distance between the cylinder and the breaking point; h b is the water depth at the breaking point; and D is the cylinder diameter. Based on the parameters used in Chella et al. [11], a numerical model was implemented. Then, the time histories of breaking wave forces, F, acting on the cylinder and the peak wave forces, F p , appearing at different cylinder positions have been calculated. Note that the forces on piles were obtained through integrating the local pressure on the wetted pile surface in time. Comparisons of the calculated results by the model presented herein versus the results of Chella et al. [11] are provided in Figure 3, where F and F p are normalized by ρgD 3 [7,8,11].
sonably satisfactory. It is only near the free surface where the calculated velocities are larger than the experimental data, reflecting the evident characteristic of wave breaking.

Comparison with Numerical Results on Breaking Wave Forces Acting on a Vertical Cylinder
Chella et al. [11] simulated breaking wave forces at different cylinder positions relative to the breaking point on an impermeable slope, and analyzed the effects of different wave shapes on the wave forces. The normalized relative distance Lc, used to distinguish the different relative positions of cylinders, is defined as: where Hb is the breaker height; xc is the distance between the cylinder and the breaking point; hb is the water depth at the breaking point; and D is the cylinder diameter. Based on the parameters used in Chella et al. [11], a numerical model was implemented. Then, the time histories of breaking wave forces, F, acting on the cylinder and the peak wave forces, Fp, appearing at different cylinder positions have been calculated. Note that the forces on piles were obtained through integrating the local pressure on the wetted pile surface in time. Comparisons of the calculated results by the model presented herein versus the results of Chella et al. [11] are provided in Figure 3, where F and Fp are normalized by ρgD 3 [7,8,11]. It can be clearly seen that the time histories of the wave forces match well with the results of Chella et al. [11], and the peak wave force is captured accurately. When the cylinder position relative to the breaking point changes, the variation trend of the calculated peak wave forces is in acceptable agreement with the results of Chella et al. [11]. Because the breaking wave forces on a pile which is located near the breaking point can change significantly when the pile position has a slight change, the major discrepancies occur at the cases of L c = −4.24 and L c = 8.48 in Figure 3. Since it is difficult to obtain the exact parameters needed to calculate the normalized relative distance, e.g., the breaker height and the water depth at the breaking point, these discrepancies are unavoidable. It can be concluded that the present model is capable of simulating the violent process of wave breaking.

Comparison with Experimental Data on Free Surface Elevations and Velocities of a Solitary Wave Propagating over a Permeable Slope or around a Permeable Structure
A series of permeable seabed experiments were conducted by Cao et al. [41], aiming to reveal the effects of seepage on the sediment initiation in the swash zone. Wave heights and fluid velocities were recorded at specific positions when a solitary wave propagated over a permeable seabed made of biochemical cotton, as given in Figure 4. "WG" and "V" denote the positions of the wave gauges and the velocimeters, respectively. Note that there is no available velocity data at the position of WG3, restricted by the minimum water depth requirement of the velocimeter.
A series of permeable seabed experiments were conducted by Cao et al. [41], aiming to reveal the effects of seepage on the sediment initiation in the swash zone. Wave heights and fluid velocities were recorded at specific positions when a solitary wave propagated over a permeable seabed made of biochemical cotton, as given in Figure 4. "WG" and "V" denote the positions of the wave gauges and the velocimeters, respectively. Note that there is no available velocity data at the position of WG3, restricted by the minimum water depth requirement of the velocimeter. A numerical wave tank similar as the experimental flume of Cao et al. [41] was configured using the present permeable seabed model. The hydraulic conductivity and the porosity of the porous media used in the experiment are 0.105 m/s and 0.4, respectively [41]. Further, the equivalent median particle diameter, d50, and other input parameters required by the numerical model can be calculated utilizing Equations (15)- (19). Excellent agreement on the wave heights and the flow velocities is found in Figure 5, indicating the present model can provide a high degree of accuracy for the permeable seabed simulations. A numerical wave tank similar as the experimental flume of Cao et al. [41] was configured using the present permeable seabed model. The hydraulic conductivity and the porosity of the porous media used in the experiment are 0.105 m/s and 0.4, respectively [41]. Further, the equivalent median particle diameter, d 50 , and other input parameters required by the numerical model can be calculated utilizing Equations (15)- (19). Excellent agreement on the wave heights and the flow velocities is found in Figure 5, indicating the present model can provide a high degree of accuracy for the permeable seabed simulations. Lara et al. [30] carried out a series of tests on the interaction between solitary waves and a porous prism on a horizontal bottom. The porous prism was located at the middle of the flume. Its d50 and porosity were 8.3 mm and 0.48, respectively. The drag coefficients can be calculated from Equations (15)- (18). A solitary wave generated at the inlet traveled towards the prism. Twelve wave gauges were installed around the prism to record the time histories of the free surface elevations. Based on the experimental setup of Lara et al. [30], a numerical wave tank was modeled.
The reliability of the permeable seabed model is the key factor to simulate the hindering effect caused by the permeable seabed accurately [42]. Thus, after waves traveling through the porous prism, the time histories of the free surface elevations around the porous prism can reflect the model's capacity for simulating the porous media. As shown in Figure 6, the calculated results are compared with the experimental data. Note that Gauge 1-2 were located before the prism and Gauge 3-6 around the prism. It can be found that the wave crest passes by the prism twice because of the wave reflection. The overall pro- Lara et al. [30] carried out a series of tests on the interaction between solitary waves and a porous prism on a horizontal bottom. The porous prism was located at the middle of the flume. Its d 50 and porosity were 8.3 mm and 0.48, respectively. The drag coefficients can be calculated from Equations (15)- (18). A solitary wave generated at the inlet traveled towards the prism. Twelve wave gauges were installed around the prism to record the time histories of the free surface elevations. Based on the experimental setup of Lara et al. [30], a numerical wave tank was modeled.
The reliability of the permeable seabed model is the key factor to simulate the hindering effect caused by the permeable seabed accurately [42]. Thus, after waves traveling through the porous prism, the time histories of the free surface elevations around the porous prism can reflect the model's capacity for simulating the porous media. As shown in Figure 6, the calculated results are compared with the experimental data. Note that Gauge 1-2 were located before the prism and Gauge 3-6 around the prism. It can be found that the wave crest passes by the prism twice because of the wave reflection. The overall process is captured precisely, and the free surface elevations around the prism are well presented in the numerical simulation. Therefore, the present model can simulate the interaction between waves and porous structures satisfactorily. of the flume. Its d50 and porosity were 8.3 mm and 0.48, respectively. The drag coefficients can be calculated from Equations (15)- (18). A solitary wave generated at the inlet traveled towards the prism. Twelve wave gauges were installed around the prism to record the time histories of the free surface elevations. Based on the experimental setup of Lara et al. [30], a numerical wave tank was modeled.
The reliability of the permeable seabed model is the key factor to simulate the hindering effect caused by the permeable seabed accurately [42]. Thus, after waves traveling through the porous prism, the time histories of the free surface elevations around the porous prism can reflect the model's capacity for simulating the porous media. As shown in Figure 6, the calculated results are compared with the experimental data. Note that Gauge 1-2 were located before the prism and Gauge 3-6 around the prism. It can be found that the wave crest passes by the prism twice because of the wave reflection. The overall process is captured precisely, and the free surface elevations around the prism are well presented in the numerical simulation. Therefore, the present model can simulate the interaction between waves and porous structures satisfactorily.  [30]. Reproduced from Lara et al. [30]; with permission from Elsevier, 2021.

Model Setup
The numerical model used in this study is illustrated in Figure 7. The origin of coordinates is located at the middle of the shoreline in the transverse direction. x, y, and z are the horizontal, transverse, and vertical coordinates, respectively. The initial water depth is set as 0.205 m and the pile diameter as 0.06 m. When the beach slope changes, the length and height of the computational domain will also be adjusted to maintain the computational efficiency of the simulations, while the width of the computational domain is always kept as 0.5 m. Note that the parameters of this numerical wave tank are similar with those in Mo et al. [7], and the results analyzed in the following sections are generated in this scaled numerical wave tank.  [30]. Reproduced from Lara et al. [30]; with permission from Elsevier, 2021.

Model Setup
The numerical model used in this study is illustrated in Figure 7. The origin of coordinates is located at the middle of the shoreline in the transverse direction. x, y, and z are the horizontal, transverse, and vertical coordinates, respectively. The initial water depth is set as 0.205 m and the pile diameter as 0.06 m. When the beach slope changes, the length and height of the computational domain will also be adjusted to maintain the computational efficiency of the simulations, while the width of the computational domain is always kept as 0.5 m. Note that the parameters of this numerical wave tank are similar with those in Mo et al. [7], and the results analyzed in the following sections are generated in this scaled numerical wave tank.
are the horizontal, transverse, and vertical coordinates, respectively. The initial water depth is set as 0.205 m and the pile diameter as 0.06 m. When the beach slope changes, the length and height of the computational domain will also be adjusted to maintain the computational efficiency of the simulations, while the width of the computational domain is always kept as 0.5 m. Note that the parameters of this numerical wave tank are similar with those in Mo et al. [7], and the results analyzed in the following sections are generated in this scaled numerical wave tank. A total of 190 cases were conducted to investigate the effect of different parameters on breaking wave forces, as listed in Table 1. Note that the water depth at the pile position hp (hp = |X·tanθ|, where X is the horizontal coordinate of the pile) is introduced to represent the pile position. All the conditions listed in Table 1 were investigated over permeable seabeds with five kinds of permeability, i.e., 9.9 × 10 −6 m 2 , 8.9 × 10 −7 m 2 , 9.9 × 10 −8 m 2 , 8.9 × 10 −9 m 2 , and 9.9 × 10 −10 m 2 . This permeability range is equivalent to the permeability of coarse sand bed or gravel bed, meaning that the compressibility of the seabed is relatively small and the seabed deformation has limited influence on the flow field. Therefore, the permeable seabed is regarded as a rigid seabed.
Solitary waves are generated on the inlet wavemaker boundary, i.e., the left boundary of the computational domain shown in Figure 7, and propagate towards the A total of 190 cases were conducted to investigate the effect of different parameters on breaking wave forces, as listed in Table 1. Note that the water depth at the pile position h p (h p = |X·tanθ|, where X is the horizontal coordinate of the pile) is introduced to represent the pile position. All the conditions listed in Table 1 were investigated over permeable seabeds with five kinds of permeability, i.e., 9.9 × 10 −6 m 2 , 8.9 × 10 −7 m 2 , 9.9 × 10 −8 m 2 , 8.9 × 10 −9 m 2 , and 9.9 × 10 −10 m 2 . This permeability range is equivalent to the permeability of coarse sand bed or gravel bed, meaning that the compressibility of the seabed is relatively small and the seabed deformation has limited influence on the flow field. Therefore, the permeable seabed is regarded as a rigid seabed. Solitary waves are generated on the inlet wavemaker boundary, i.e., the left boundary of the computational domain shown in Figure 7, and propagate towards the positive direction of the x axis. The Sommerfeld radiation condition applied on the outlet boundary ( Figure 7) allows fluid to flow away from the computational domain without significant reflection. That helps to reduce the length of the computational domain and improve the computational efficiency. A no-slip condition is imposed to treat the bottom boundary and the pile surface. Other boundaries are symmetry conditions [10].

Mesh Dependency Analysis
The effect of mesh dependency on the accuracy of the numerical simulations is of vital importance. To investigate the mesh dependency in the present model, a grid convergence study whose conditions are similar with the experimental setup of Mo et al. [7] was conducted. The computational domain is 10 m long, 0.5 m wide, and 0.8 m high. The slope angle, pile diameter, initial water depth and initial wave height are 5.1 • , 0.06 m, 0.205 m, and 0.067 m, respectively, also as shown in Figure 7.
The computational domain is partitioned into structural meshes with regular-hexahedron shapes, and the grid sizes are uniform in all directions. Six different grid sizes adopted in the simulations, and the corresponding total grid number in the computational domain, are shown in Table 2. The comparison of free surface elevation time series and the maximum run-up height for different grid sizes is given in Figure 8, where all the heights are normalized by h 0 . Note that in Figure 8b the dotted line corresponds to the experimental result in Mo et al. [7].   Figure 9 shows the variation of the normalized peak wave forces on a pile with a solitary wave propagating with different slope angles on a permeable seabed in the representative case (H0 = 0.067 m, hp = 23.20 mm, α = 0°). It can be found that the variation trends of peak wave forces are quite similar on different permeable seabeds, i.e., firstly increasing and then decreasing until stabilization. The peak wave force reaches the maximum for the 5.1° slope with the permeability of K = 9.9 × 10 −10 m 2 , because the wave breaks right before the pile at this slope, as shown in Figure 10d. That coincides with the result of previous studies conducted on impermeable slopes [11]. It's clear that the changes of the grid size have minimum effects on the wave height when solitary waves are propagating on a plane seabed. The calculated wave shape coincides well with the experimental data. However, after a wave propagates along the slope and slams against the pile, the maximum run-up height calculated with a coarse grid shows significant discrepancies with the experimental data. The numerical results are verging on stabilization when the grid size is smaller than 5 mm, as shown in Figure 8. Consequently, the grid size of dx = 5 mm is adopted in all simulations in terms of the accuracy of numerical results and the constrains of computational capacity. Figure 9 shows the variation of the normalized peak wave forces on a pile with a solitary wave propagating with different slope angles on a permeable seabed in the representative case (H 0 = 0.067 m, h p = 23.20 mm, α = 0 • ). It can be found that the variation trends of peak wave forces are quite similar on different permeable seabeds, i.e., firstly increasing and then decreasing until stabilization. The peak wave force reaches the maximum for the 5.1 • slope with the permeability of K = 9.9 × 10 −10 m 2 , because the wave breaks right before the pile at this slope, as shown in Figure 10d. That coincides with the result of previous studies conducted on impermeable slopes [11].

Slope Angles
solitary wave propagating with different slope angles on a permeable seabed in the representative case (H0 = 0.067 m, hp = 23.20 mm, α = 0°). It can be found that the variation trends of peak wave forces are quite similar on different permeable seabeds, i.e., firstly increasing and then decreasing until stabilization. The peak wave force reaches the maximum for the 5.1° slope with the permeability of K = 9.9 × 10 −10 m 2 , because the wave breaks right before the pile at this slope, as shown in Figure 10d. That coincides with the result of previous studies conducted on impermeable slopes [11].   Figure 10 also shows that on the slope with an angle of θ = 3°, the wave breaks much before interacting with the pile. And the type of breaker is more likely to be spilling instead of plunging, leading to a remarkable decrease in breaking wave forces. On the slopes Figure 10. Velocity fields around the pile when the peak wave force on the pile appears on a permeable seabed with different slope angles (unit: m/s). Figure 10 also shows that on the slope with an angle of θ = 3 • , the wave breaks much before interacting with the pile. And the type of breaker is more likely to be spilling instead of plunging, leading to a remarkable decrease in breaking wave forces. On the slopes whose angles are larger than 5.1 • , waves don't break before the pile and keep similar shapes until slamming against the pile. Thus, the peak wave force becomes relatively small and shows little difference as the slope angles increase.
With the seabed permeability increasing, the breaking point shifts shoreward due to the wave absorbing effect of a permeable seabed. That results in the decrease in the maximum peak wave force on the 5.1 • slope. Whereas, the peak wave force on the 5.1 • slope is still the maximum even when the slope permeability is the largest (K = 9.9 × 10 −6 m 2 ). It turns out that slope angles play a dominant role in influencing the maximum breaking wave force compared with the seabed permeability in present cases.

Wave Heights
Since wave heights could influence the breaker type directly, it is evident that wave force generally increases as the relative wave height (H 0 /h) increases [43].
The normalized peak wave forces under different wave heights, 0.03 m and 0.09 m, respectively, are shown in Figure 11. The columns denote the magnitude of the increase in the peak wave force when the wave height increases from 0.03 m to 0.09 m. Obviously, there is a growing trend in the magnitude with seabed permeability decreasing, since less energy is dissipated by the viscous friction inside the permeable seabed. But an exception is observed on the slope with an angle of θ = 5.1 • . in the peak wave force when the wave height increases from 0.03 m to 0.09 m. Obviously, there is a growing trend in the magnitude with seabed permeability decreasing, since less energy is dissipated by the viscous friction inside the permeable seabed. But an exception is observed on the slope with an angle of θ = 5.1°. Figure 11. The magnitude of the increase in the peak wave forces caused by the increase of wave heights on the permeable seabed with different slope angles. Figure 12 illustrates that waves with small heights (H = 0.03 m) do not break before the pile, which makes the Fp (represented by triangles in Figure 11) remain relatively stable. Thus, the variation trend of ∆Fp mainly depends on the Fp which is represented by rectangles in Figure 11. On the 5.1° slope, the wave breaks slightly before the pile with a wave height of H = 0.09 m and seabed permeability of K = 9.9 × 10 −6 m 2 (Figure 12b), leading to a comparatively large peak wave force.
However, as seabed permeability decreases to 9.9 × 10 −8 m 2 , waves propagating over the permeable seabed have already broken before slamming against the pile. The overturning wave crest interacting with the pile resulted in significant decrease in the peak wave forces. While seabed permeability decreases further to 9.9 × 10 −10 m 2 , the broken wave crest hit the still water before the pile and cause a more significant energy loss around the wave crest. But it is clearly visible that the broken crest also produces a "secondary wave wall" upstream of the pile, as shown in Figure 12f. The secondary wave wall slams against piles immediately following the wave crest, generating a larger peak wave force even than on the slope with permeability of K = 9.9 × 10 −6 m 2 .
In addition, with the slope angles increasing, the magnitude of the increase caused by wave breaking decreases. As proved in Section 3.1, it can be attributed to waves interacting with pile earlier (Figure 12g,h). Figure 11. The magnitude of the increase in the peak wave forces caused by the increase of wave heights on the permeable seabed with different slope angles. Figure 12 illustrates that waves with small heights (H = 0.03 m) do not break before the pile, which makes the F p (represented by triangles in Figure 11) remain relatively stable. Thus, the variation trend of ∆F p mainly depends on the F p which is represented by rectangles in Figure 11. On the 5.1 • slope, the wave breaks slightly before the pile with a wave height of H = 0.09 m and seabed permeability of K = 9.9 × 10 −6 m 2 (Figure 12b), leading to a comparatively large peak wave force.
However, as seabed permeability decreases to 9.9 × 10 −8 m 2 , waves propagating over the permeable seabed have already broken before slamming against the pile. The overturning wave crest interacting with the pile resulted in significant decrease in the peak wave forces. While seabed permeability decreases further to 9.9 × 10 −10 m 2 , the broken wave crest hit the still water before the pile and cause a more significant energy loss around the wave crest. But it is clearly visible that the broken crest also produces a "secondary wave wall" upstream of the pile, as shown in Figure 12f. The secondary wave wall slams against piles immediately following the wave crest, generating a larger peak wave force even than on the slope with permeability of K = 9.9 × 10 −6 m 2 .
around the wave crest. But it is clearly visible that the broken crest also produces a "secondary wave wall" upstream of the pile, as shown in Figure 12f. The secondary wave wall slams against piles immediately following the wave crest, generating a larger peak wave force even than on the slope with permeability of K = 9.9 × 10 −6 m 2 .
In addition, with the slope angles increasing, the magnitude of the increase caused by wave breaking decreases. As proved in Section 3.1, it can be attributed to waves interacting with pile earlier (Figure 12g,h).

Pile Positions and Inclination Angles
As has been noted, the process of wave breaking could be influenced by slope angles and incident wave heights. In order to investigate the effects of the pile position on the breaking wave force, the pile position L, the distance between the pile and the shoreline, are normalized by (H0/tanθ) to exclude the influence of other parameters as follows: where X0 is the  In addition, with the slope angles increasing, the magnitude of the increase caused by wave breaking decreases. As proved in Section 3.1, it can be attributed to waves interacting with pile earlier (Figure 12g,h).

Pile Positions and Inclination Angles
As has been noted, the process of wave breaking could be influenced by slope angles and incident wave heights. In order to investigate the effects of the pile position on the breaking wave force, the pile position L, the distance between the pile and the shoreline, are normalized by (H 0 /tanθ) to exclude the influence of other parameters as follows: where X 0 is the  [8,11]. Hence these cases are not investigated here. The effect of the pile positions relative to the breaking point on the breaking wave force has been preliminarily investigated by Liu et al. [44]. It has been found that the positions where the pile suffers maximum wave force can transfer from the breaking point to the position 1D behind the breaking point, when the seabed permeability increases from 9.9 × 10 −10 m 2 to 8.9 × 10 −7 m 2 . On the high permeability seabed (K = 9.9 × 10 −6 m 2 ), breaking wave forces on the pile have limited changes at various pile positions, due to the enhanced wave damping effect of the permeable seabed. However, the changes of the pile inclination angles have not been considered in Liu et al. [44] Figure 13 shows the normalized peak wave forces with different pile positions and inclined angles founded on a permeable seabed. When the pile is relatively far from the shoreline (L = −0.70), the breaking wave force on the pile has limited changes as the inclined angle of the pile differs from −45 • to 45 • . That is because the wave hasn't started breaking and the wave profile remains nearly stable before the wave hit the pile. to the position 1D behind the breaking point, when the seabed permeability increases from 9.9 × 10 −10 m 2 to 8.9 × 10 −7 m 2 . On the high permeability seabed (K = 9.9 × 10 −6 m 2 ), breaking wave forces on the pile have limited changes at various pile positions, due to the enhanced wave damping effect of the permeable seabed. However, the changes of the pile inclination angles have not been considered in Liu et al. [44] Figure 13 shows the normalized peak wave forces with different pile positions and inclined angles founded on a permeable seabed. When the pile is relatively far from the shoreline (L′ = −0.70), the breaking wave force on the pile has limited changes as the inclined angle of the pile differs from −45° to 45°. That is because the wave hasn't started breaking and the wave profile remains nearly stable before the wave hit the pile.  For piles installed around the wave breaking point (L = −0.43-−0.27), the peak wave force generally increases firstly and then decreases with the inclined angles increasing. Similar with L = −0.70, the pile with an inclined angle of α = −45 • doesn't suffer large breaking wave forces. The maximum peak wave force occurs at α = −22.5 • or α = 0 • , namely when waves break right on the front surface of the pile, as Figure 14a-d show. Sometimes, the waves that have already broken slightly could also produce large breaking wave forces on the pile with an inclination angle of α = 45 • (Figure 14e

Conclusions
In this study, a numerical wave flume which contained a slope seabed was simulated. The Forchheimer saturated drag equation was introduced to consider the permeability of the seabed and simulate the effect of porous media on fluid motions. Verifications with previous experimental and numerical data showed good agreement. Then, 190 numerical simulations were conducted and the characteristics of breaking solitary wave forces on piles founded on a permeable seabed were analyzed systematically. The main conclusions can be summarized as follows: 1. The slope angles play a dominant role in influencing breaking wave forces on a pile founded on a permeable seabed. The maximum breaking wave force appears when waves break just before the pile. Once waves have broken, the larger slope angles will have a reduced effect on the breaking wave force. 2. Breaking wave forces increase with the incident wave height increasing. However, Figure 14. Velocity fields around the pile when the peak wave force on the pile appears on a permeable seabed with different pile positions and inclined angles (unit: m/s).

Conclusions
In this study, a numerical wave flume which contained a slope seabed was simulated. The Forchheimer saturated drag equation was introduced to consider the permeability of the seabed and simulate the effect of porous media on fluid motions. Verifications with previous experimental and numerical data showed good agreement. Then, 190 numerical simulations were conducted and the characteristics of breaking solitary wave forces on piles founded on a permeable seabed were analyzed systematically. The main conclusions can be summarized as follows: 1.
The slope angles play a dominant role in influencing breaking wave forces on a pile founded on a permeable seabed. The maximum breaking wave force appears when waves break just before the pile. Once waves have broken, the larger slope angles will have a reduced effect on the breaking wave force.

2.
Breaking wave forces increase with the incident wave height increasing. However, the magnitude of the increase generally decreases as the slope angles and the permeability increase.

3.
On the permeable seabed, extreme breaking wave forces can also occur when a "secondary wave wall" interacts with the pile after wave breaking.

4.
Breaking wave forces acting on a pile can be influenced by the inclination angle of the pile significantly, especially when the pile was installed around the wave breaking point. The maximum peak wave forces usually occur at α = −22.5 • or α = 0 • . Results of the variation trends of breaking wave force with different parameters related to waves, permeable seabed and piles, can be used to extend the knowledge in designing the pile foundation or the high-rise pile cap foundation of OWTs and assessing the potential risk of coastal infrastructure.
Author Contributions: Z.L., methodology, investigation, writing-original draft preparation; Z.G., conceptualization, methodology, supervision, writing-review and editing; Y.D., visualization; F.Z., visualization. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request.