Boussinesq Simulation of Coastal Wave Interaction with Bottom-Mounted Porous Structures

: A Boussinesq-type wave model is developed in this paper to simulate the interaction of coastal waves with bottom-mounted porous structures. The governing equations are rewritten in the conservative form to facilitate the use of hybrid ﬁnite volume (FV) and ﬁnite difference (FD) method. Higher-order slope terms are also inserted into the equations to account for rapidly varying bathymetry. The convective ﬂux is approximated using the FV method, while the remaining terms are discretized using the FD method in a uniform rectangle grid system. The time integration is implemented using the third order Runge–Kutta method with an adaptive time step. A single GPU parallel computation is also implemented to save computation costs. The numerical model is validated against a series of experimental datasets, including data acquired in a new laboratory experiment. The predictions are in overall agreement with the measurements, proving that the model is capable of handling wave interaction with porous structures in the coastal region for a wide range of scenarios.


Introduction
Artificial porous structures such as rubble-mound breakwaters and revetments are commonly built to protect harbors, coastal buildings, beaches and shorelines from erosion. These kinds of coastal structures partially reflect and partially dissipate wave energy through wave-structures drag resistance and inertial resistance force, resulting in only part of the wave transmission. The interaction of nonlinear, shallow water waves with porous coastal structures is therefore an important subject in coastal planning and design.
With the aid of rapid advancements in computing power, numerical modeling approaches have been becoming increasingly popular when it comes to studying the hydrodynamics involved in wave interaction with porous structures. The state-of-the-art models are those based on the Reynolds averaged Navier-Stokes (RANS) equations (e.g., [1][2][3][4] and references therein). RANS-type models employ few simplifying assumptions compared to other depth-integrated models, such as nonlinear shallow water equations, e.g., [5,6] and Boussinesq-type equations, e.g., [7][8][9][10]. As a result, these models can predict the details of wave-structure interaction processes. However, for three-dimensional phenomena such as wave refraction and diffraction, solving the Navier-Stokes equations takes a long time. Thus, these models are rarely used in practical applications. In recent years, a large amount of progress has also been achieved using Smoothed Particle Hydrodynamics (SPH) models, which are meshless and have comparable behaviors to RANS-type models (e.g., [11][12][13]). These models still suffer from a huge computational burden that limits their practical applications for modeling wave interaction with porous structures.

Governing Equations
Taking inertial and drag resistances into account, Vu et al. [10] have recently derived a set of Boussinesq equations for wave propagation in porous medium, which reduces to the equations of Madsen and Sørensen [25] in a non-porous medium. In their work, the equations are in non-conservative forms and are solved numerically using the FD scheme. To facilitate the use of the FV scheme, the Boussinesq equations are rewritten in the conservative form after some mathematical manipulations (see [25,30]) under the Cartesian coordinates ( Figure 1). ∂η ∂t with ψ x = β −(B + 1 3 )h 2 P xxt + Q xyt − hh x ( 1 3 P xt + 1 6 Q yt ) −hh y 1 6 Q xt − 1 6 hh xx P t − 1 6 hh xy Q t + 1 3 h 2 x P t + 1 3 h x h y Q t −Bgh 3 (η xxx + η xyy ) −hh x (2Bghη xx + Bghη yy ) − hh y (Bghη xy ) −hh xx Bghη x − hh xy Bghη y +α P − (B + 1 3 )h 2 P xx + Q xy − hh x ( 1 3 P x + 1 6 Q y ) − hh y 1 6 Q x − 1 6 hh xx P − 1 6 hh xy Q + 1 3 h 2 x P + 1 3 h x h y Q ψ y = β −(B + 1 3 )h 2 Q yyt + P xyt − hh y ( 1 3 Q yt + 1 6 P xt ) −hh x 1 6 P yt − 1 6 hh yy Q t − 1 6 hh xy P t + 1 3 h 2 y Q t + 1 3 h x h y P t −Bgh 3 (η yyy + η xxy ) −hh y (2Bghη yy + Bghη xx ) − hh x (Bghη xy ) −hh yy Bghη y − hh xy Bghη x +α Q − (B + 1 3 )h 2 Q yy + P xy − hh y ( 1 3 Q y + 1 6 P x ) − hh x 1 6 P y − 1 6 hh yy Q − 1 6 hh xy P + 1 where η is the free surface elevation measured from the still water level, h is the still water depth, d = h + η is the total water depth, g is the gravity acceleration, P = du and Q = dv are fluxes with u = (u, v) is the depth-averaged velocity vector and B = 1/15 is the dispersion parameter. β = λ + (1 − λ)(1 + κ) is the inertial resistance coefficient, with λ and κ = 0.34 being the porosity of porous medium and the added mass coefficient; α is the drag resistance coefficient calculated by where ν is the kinematic viscosity of water, and d 50 is the size of the solid material. α l and α t are coefficients representing the laminar and turbulent flow resistances, respectively. The values of these two parameters depend on the porous material property, the Reynolds number and the flow characteristics. The universal values are not available, and we calibrate them using the experimental data (provided in Section 4) following the overwhelming majority of previous studies, e.g., [1][2][3][4][5][6][7][8][9][10].   It should be mentioned here that (i) λ = 1 and λ < 1 in Equations (1)-(3) represent wave motion in porous and non-porous mediums, respectively, which simplifies the numerical scheme, as noted before. Otherwise, additional efforts should be made to deal with the matching conditions along the interface between porous and non-porous mediums (see [9,23]); (ii) Higher-order slope terms, originally derived by Kim et al. [30], are inserted into Equations (1)-(3) to improve model performance for rapidly varying topography.

Numerical Scheme and Boundary Conditions
The hybrid FV/FD shock-capturing scheme, originally proposed by the authors for wave propagation in open water [27], is adopted to solve Equations (1)-(3). These equations are almost identical to those presented in [27], except for extra terms related to porous medium, and the numerical implementation is primarily straightforward. Furthermore, the methods employed in [28] are used to deal with moving shoreline and breaking waves, which are not well addressed in [27]. The main procedures are summarized herein for clarity and completeness.
On a rectangular cell system, the convective flux terms are solved by the FV method, and the rest of the terms in the equations are discretized using the FD method. The hydrostatic construction technique is used during the implementation of the fourth-order MUSCL (Monotone Upstream-Centred Schemes for Conservation Laws) with the Minmod limiter function, which helps to capture the moving shoreline in an accurate and efficient manner [28]. A central upwind scheme is then used to calculate the interface fluxes. The third-order Runge-Kutta scheme with the adaptive time step is applied to performing time integration. The time step ∆t is restricted by the Courant number C r , which is set to 0.20 for all simulations. The velocity components are obtained by solving the tridiagonal system, resulting from discretizing the time derivative terms in Equations (2) and (3) with the second-order finite difference formula [27]. The specific values of grid size ∆x will be mentioned explicitly in each test case, which are found to be sufficient to yield numerical stability and convergence for all the simulations in this study.
The unique feature of a shock-capturing Boussinesq model is its ability to capture breaking waves as discontinuities by deactivating dispersive terms. In the computation, the ratio of wave height to water depth ε (estimated approximately by η/h) and the local wave angle ∇η are computed for each cell. If at least one of the criteria is satisfied, the corresponding cells are labeled as breaking waves and dispersive terms are deactivated. Dispersive terms can only be reactivated once ε reduces to below a certain threshold ε < 0.35. This hybrid wave breaking criteria does not need to track breaking fronts and admits the steady jumps or initial discontinuities [28].
Three ghost cells are employed to enclose the entire computational domain. For the reflective wall, the tangential and normal velocity components on ghost cells are determined from the inner domain by imposing symmetric and anti-symmetric conditions with respect to the solid wall [27]. For incident wave boundary, the variables on the ghost cells are given according to the theoretical wave solutions, i.e., solitary wave theory, linear wave theory and cnoidal wave theory in the present study, as known quantities. In this way, wave signals are sent into the computational domain to mimic wave generation from a paddle-type wavemaker. A spongey layer is also placed at the front of solid walls to absorb wave energy once necessary [14].

GPU Implementation
Boussinesq equations are phase resolving, and usually require a large number of discrete grids to obtain detailed intra-wave hydrodynamics and high resolution of coastal structures. To meet the stability requirement, the time step is usually in the order of O(T/100, where T is the typical wave period). Thus, BT simulation is computationally demanding at fine resolutions in practical applications, although it is much cheaper than RANS simulation. Recently, Fang et al. [29] successfully ported their BT model [27] to a single GPU via CUDA C. The massive computation tasks (e.g., fluxes evaluation, MUSCL reconstruction) were coded as kernel functions and executed on the GPU device. To achieve better performance, the model uses a CR (Cyclic Reduction) technique to solve massive tridiagonal linear systems and overlapped tiling/shared memory to reduce global memory access and enhance data reuse. Compared to the serial code, the maximum speedup ratios for single-and double-precision calculations are 55.56 and 32.57, respectively for tested cases. We adopt the same procedure for GPU parallel computation, further details of which can be found in [29]. High computation efficiency is another advantage of the present BT model, and to the authors' knowledge, only a few BT models have been accelerated by the GPU technique to date [31][32][33].

Physical Wave Flume Experiments
Previous laboratory datasets used to validate the BT model mainly concern solitary wave and statistical wave height of short waves, while time series of the free surface elevation of dispersive waves are less reported. New laboratory experiments are thus conducted to provide a larger dataset. Figure

Dam Break on a Porous Structure
The first test case is to simulate the experiments of dam break on a porous structure conducted by Liu et al. [1]. The water tank used in the experiments is 0.892 m long, 0.44 m wide, and 0.58 m high. Crushed rocks with a mean diameter of d50 = 0.0159 m and porosity λ = 0.49 were placed in the center of the water tank to form a 0.29 m long and 0.58 m tall dam. The reservoir has an initial water height of 0.25 m and is separated from the porous dam by a 0.02 m thick gate. The numerical setup is the same as the experimental setup and is illustrated in the first panel of Figure 3, which corresponds to the initial state at t = 0 s. The simulation domain is discretized into cells of Δx = 0.005 m. The drag coefficients of αl = 4000, αt = 4.2 are calibrated by running trial-and-error simulations to achieve the best agreement between numerical results and experimental data.
To evaluate the agreements between numerical results and experimental data for a given quantity v, the Willmott index [34] is used. This index is introduced as where x(j) is the measured data point, y(j) is the computation result, and x is the mean value of series y(j). Perfect agreement is indicated by W = 1, whereas W = 0 indicates complete disagreement. Figure 3 compares the simulated water level with the experimental data. The overall processes such as upward jet as the water body crashes into the porous medium (t = 0.2

Dam Break on a Porous Structure
The first test case is to simulate the experiments of dam break on a porous structure conducted by Liu et al. [1]. The water tank used in the experiments is 0.892 m long, 0.44 m wide, and 0.58 m high. Crushed rocks with a mean diameter of d 50 = 0.0159 m and porosity λ = 0.49 were placed in the center of the water tank to form a 0.29 m long and 0.58 m tall dam. The reservoir has an initial water height of 0.25 m and is separated from the porous dam by a 0.02 m thick gate. The numerical setup is the same as the experimental setup and is illustrated in the first panel of Figure 3, which corresponds to the initial state at t = 0 s. The simulation domain is discretized into cells of ∆x = 0.005 m. The drag coefficients of α l = 4000, α t = 4.2 are calibrated by running trial-and-error simulations to achieve the best agreement between numerical results and experimental data.

Solitary Wave Passing through Porous Breakwaters
Two cases from the new experiments (S1 and S2 in Table 1) are first simulated to show the model performance in terms of solitary wave interaction with porous structure. The numerical flume is the same as the physical flume shown in Figure 2. The grid size is 0.01687 m, and αl = 800, αt = 9.0. A 3.0 m-wide sponger layer is placed at the end of the numerical flume to absorb wave energy. Figures 4 and 5 show the comparison of the computed free surface elevation with the experiment data, respectively. For case S1, the wave gauges located in front of the porous structure (G1-G3) clearly record two peaks, corresponding to the incident and reflected waves from the porous structure. The BT model captures the arrival time of both reflected and transmitted waves, but it overestimates the amplitude of the secondary wave at G1 and G2. Adjacent to the porous wall at G5, these two peaks merge into one with a large amplitude, about 1.5 times that of the incident amplitude. The solitary wave shape is well maintained on the leeside of the breakwater at G6-G8, where the wave amplitude is largely dissipated by 71% due to the presence of porous breakwater. Case S2 is more non- To evaluate the agreements between numerical results and experimental data for a given quantity v, the Willmott index [34] is used. This index is introduced as where x(j) is the measured data point, y(j) is the computation result, and x is the mean value of series y(j). Perfect agreement is indicated by W = 1, whereas W = 0 indicates complete disagreement. Figure 3 compares the simulated water level with the experimental data. The overall processes such as upward jet as the water body crashes into the porous medium (t = 0.2 s), the reflection of dam break flow from the right end of water tank (t = 1.4~1.8 s) and the trend of water seeping porous structure are reasonably predicted by the BT model. It is noted that water close to the bottom moves earlier than the water at the free surface when the gate is opened manually in the experiments. However, the depth-integrated BT model cannot reproduce this process and leads to significant inconsistencies with the measurements at t = 0.2, 0.4, 0.6 and 0.8 s. Better agreements are achieved as the flow evolves (t = 1.0-2.2 s). The RANS simulation results using REEF3D in [4] are also shown in Figure 3 as a comparison. The mean value of the Willmott index from the present simulation and RANS simulation are 0.953 and 0.962, respectively. RANS simulations present better predictions, as expected, except for the disagreements at the initial stage due to the reasons mentioned above. After t = 0.8 s, the simulation results of the present BT model are comparable to the predictions of RANS simulations. BT models based on the FD method fail to simulate this test case, since they are unable to treat the water discontinuity at t = 0 s, which demonstrates the merits of the shock-capturing hybrid scheme adopted in the present study.

Solitary Wave Passing through Porous Breakwaters
Two cases from the new experiments (S1 and S2 in Table 1) are first simulated to show the model performance in terms of solitary wave interaction with porous structure. The numerical flume is the same as the physical flume shown in Figure 2. The grid size is 0.01687 m, and α l = 800, α t = 9.0. A 3.0 m-wide sponger layer is placed at the end of the numerical flume to absorb wave energy. Figures 4 and 5 show the comparison of the computed free surface elevation with the experiment data, respectively. For case S1, the wave gauges located in front of the porous structure (G1-G3) clearly record two peaks, corresponding to the incident and reflected waves from the porous structure. The BT model captures the arrival time of both reflected and transmitted waves, but it overestimates the amplitude of the secondary wave at G1 and G2. Adjacent to the porous wall at G5, these two peaks merge into one with a large amplitude, about 1.5 times that of the incident amplitude. The solitary wave shape is well maintained on the leeside of the breakwater at G6-G8, where the wave amplitude is largely dissipated by 71% due to the presence of porous breakwater. Case S2 is more nonlinear than Case S1, but a similar evolution process is observed. For this case, the dimensionless transmission amplitude is reduced to 0.21, which is inconsistent with previous findings that large nonlinearity causes a smaller transmission coefficient [23,35]. In summary, the waveform and phase are predicted very well for both cases.
Regarding a solitary wave passing through porous structures, previous studies in the literature mainly concerned with the transmitted and reflected wave amplitude. We further simulate six cases from the experiments of Vidal et al. [35], where transmission and reflection coefficients are available and have been extensively used for model validation. Table 2 summarizes the material property and dimension of porous breakwaters from the experiments. In the simulation, a grid size of 0.025 m is used, and the drag coefficients are set to α l = 800, α t = 3.0. Figure 6 compares the transmission coefficients K t from the current simulation with the experimental data, and they are in good agreement. The reflection coefficients K r are also plotted in the figure. As nonlinearity increases, the decreasing trend of K t and the increasing trend of K r are well predicted by the BT simulation. linear than Case S1, but a similar evolution process is observed. For this case, the dimensionless transmission amplitude is reduced to 0.21, which is inconsistent with previous findings that large nonlinearity causes a smaller transmission coefficient [23,35]. In summary, the waveform and phase are predicted very well for both cases.

Regular Wave Passing through Porous Breakwaters
In this section, numerical simulations are conducted for regular wave cases (R1-R in Table 1) to test the model's ability to describe short dispersive wave interaction with porous structures. A grid size of Δx = 0.02 m is used to discretize the computational do main. A 3.0 m-wide sponger layer is placed at the end of the numerical flume to absorb wave energy. Figures 7-9 compare the time series of computed and measured free surface elevation at eight wave gauge locations, five in front of the porous structure and three leeward o the breakwater. For all cases, both the wave shape and phase are reproduced by the BT model with a high degree of accuracy. The incident waves are partially reflected from th porous structure, which interacts with the incident waves and creates a standing wav pattern in front of the structure (G1-G5). These wave profiles are very well caught during the simulation, including the initial state and the subsequent quasi-steady wave pattern On the leeward side of the porous breakwater, the wave profiles calculated at location G6-G8 correspond to the incident wave that is damped by the porous structure. The wav shape and phase are also predicted very well by the BT model for three cases.

Regular Wave Passing through Porous Breakwaters
In this section, numerical simulations are conducted for regular wave cases (R1-R3 in Table 1) to test the model's ability to describe short dispersive wave interaction with porous structures. A grid size of ∆x = 0.02 m is used to discretize the computational domain. A 3.0 m-wide sponger layer is placed at the end of the numerical flume to absorb wave energy. Figures 7-9 compare the time series of computed and measured free surface elevation at eight wave gauge locations, five in front of the porous structure and three leeward of the breakwater. For all cases, both the wave shape and phase are reproduced by the BT model with a high degree of accuracy. The incident waves are partially reflected from the porous structure, which interacts with the incident waves and creates a standing wave pattern in front of the structure (G1-G5). These wave profiles are very well caught during the simulation, including the initial state and the subsequent quasi-steady wave pattern. On the leeward side of the porous breakwater, the wave profiles calculated at locations G6-G8 correspond to the incident wave that is damped by the porous structure. The wave shape and phase are also predicted very well by the BT model for three cases.

Wave Passing through Porous Breakwater in a Wave Basin
In order to examine the model performance for three-dimensional problems, numerical results are compared with the experimental measurements conducted by Lara et al. [36] in a wave basin. The sketch of the experiments and 15 wave gauges are shown in Figure 10. Table 3

Wave Passing through Porous Breakwater in a Wave Basin
In order to examine the model performance for three-dimensional problems, numerical results are compared with the experimental measurements conducted by Lara et al. [36] in a wave basin. The sketch of the experiments and 15 wave gauges are shown in Figure 10. Table 3  The simulated interaction processes between solitary wave and porous structure are presented in Figure 11. As can be seen from this figure, wave damping and wave transmission resulting from wave penetration through the porous structure, wave reflection from porous structure and side walls, wave diffraction and wave run-up on structure result in a complicated wave field. The solitary wave is approaching the porous structure at t = 5.0 s and t = 6.0 s. The incident wavefront is separated by the structure and the transmitted wave in the unblocked region, as seen at t = 7.0 s. The solitary wave climbs up the structure at the same time, and the subsequent run down is observed at t = 7.25 s. The panels corresponding to t = 7.85 s and t = 8.45 s show the wave transmitted through and beyond the structure approaching the end of the wave basin and the first partially reflected wave from the structure traveling toward the wavemaker. The end wall of the wave basin reflects the waves, and these waves propagate towards the structure at t = 10.0 s and t = 11.5 s. The second time wave interaction will occur subsequently (not shown in Figure 11).   The simulated interaction processes between solitary wave and porous structure are presented in Figure 11. As can be seen from this figure, wave damping and wave transmission resulting from wave penetration through the porous structure, wave reflection from porous structure and side walls, wave diffraction and wave run-up on structure result in a complicated wave field. The solitary wave is approaching the porous structure at t = 5.0 s and t = 6.0 s. The incident wavefront is separated by the structure and the transmitted wave in the unblocked region, as seen at t = 7.0 s. The solitary wave climbs up the structure at the same time, and the subsequent run down is observed at t = 7.25 s. The panels corresponding to t = 7.85 s and t = 8.45 s show the wave transmitted through and beyond the structure approaching the end of the wave basin and the first partially reflected wave from the structure traveling toward the wavemaker. The end wall of the wave basin reflects the waves, and these waves propagate towards the structure at t = 10.0 s and t = 11.5 s. The second time wave interaction will occur subsequently (not shown in Figure 11).  A further comparison of the time series of free surface elevation at 15 wave gauges with experimental data is presented along with the RANS type simulation using REEF3D [4] in Figure 12. The incident wave and the reflected wave from the end of the wave basin can be seen at G1. G2 records the incident wave, the partially reflected wave from the porous structure and the reflected wave from the end of the wave basin in sequence. Both G3 and G4 record a secondary crest, as they are positioned just 1.0 m from the porous structure and are quickly affected by wave reflection. G5 and G6 are located around the head of the porous structure and have a similar variation trend. Wave gauges G7-G11 and G13 are placed on the leeside of the porous structure. The amplitude of the transmitted wave at these locations is greatly damped by the porous structure. In contrast, the second peak formed by the reflected wave from the end of the domain has a relatively larger amplitude. G12 is located further behind the structure head and the incident wave consists of the unaffected part of the wave and the part transmitted through the structure. It also records the reflected waves from the basin wall. G14 and G15 are fully exposed to the incident wave and are less affected by the porous structure. The current BT simulations agree well with the experimental data, with the deviations mainly observed in the phase of the reflected waves. This lag is attributed to discrepancies in several measurements, including the breakwater geometry, its installation location in the basin and slight variations in the wave gauge locations [36]. The developed BT model shows a high degree of accuracy in predicting the free surface time series with the mean value of W = 0.949, achieving almost identical performance to that of REEF3D (W = 0.952). A further comparison of the time series of free surface elevation at 15 wave gauges with experimental data is presented along with the RANS type simulation using REEF3D [4] in Figure 12. The incident wave and the reflected wave from the end of the wave basin can be seen at G1. G2 records the incident wave, the partially reflected wave from the porous structure and the reflected wave from the end of the wave basin in sequence. Both G3 and G4 record a secondary crest, as they are positioned just 1.0 m from the porous structure and are quickly affected by wave reflection. G5 and G6 are located around the head of the porous structure and have a similar variation trend. Wave gauges G7-G11 and G13 are placed on the leeside of the porous structure. The amplitude of the transmitted wave at these locations is greatly damped by the porous structure. In contrast, the second peak formed by the reflected wave from the end of the domain has a relatively larger amplitude. G12 is located further behind the structure head and the incident wave consists of the unaffected part of the wave and the part transmitted through the structure. It also The simulated three-dimensional wave field for the cnoidal wave case is plotted in Figure 13. Figure 14 compares the time series of simulated free surface elevation with the laboratory measurements and the OpenFOAM simulations [3]. Unlike the previous isolated solitary wave case, cnoidal waves are sent into the wave basin continuously. Wave reflections from porous structures and three sidewalls, the interaction of incident and reflected waves and the interaction between these waves and porous structure are thus more complicated. Overall, the BT model captures the wave shape and amplitude at all wave gauge locations with reasonable accuracy, although less satisfactory results are observed towards the end of the simulation. The capability of the BT model to predict three-dimensional dispersive problems is proved. The mean value of the Willmott index from the present simulation and RANS simulation are 0.902 and 0.906, respectively, for this case, illustrating almost equivalent performance. It is reported that 20 s are simulated in less than 17 h using 128 processors (2.6 GHz) for OpenFOAM simulations [3]. It only takes 5.3 min using two processors (3.6 GHz) with the aid of NVIDIA GeForce GTX 970 card for the present BT simulation. The simulated three-dimensional wave field for the cnoidal wave case is plotted in Figure 13. Figure 14 compares the time series of simulated free surface elevation with the laboratory measurements and the OpenFOAM simulations [3]. Unlike the previous sional dispersive problems is proved. The mean value of the Willmott index from the present simulation and RANS simulation are 0.902 and 0.906, respectively, for this case, illustrating almost equivalent performance. It is reported that 20 s are simulated in less than 17 h using 128 processors (2.6 GHz) for OpenFOAM simulations [3]. It only takes 5.3 min using two processors (3.6 GHz) with the aid of NVIDIA GeForce GTX 970 card for the present BT simulation.

Conclusions
A set of extended Boussinesq equations [10] is numerically solved using a hybrid FV/FD scheme proposed by the authors [27,29] to simulate wave propagation through

Conclusions
A set of extended Boussinesq equations [10] is numerically solved using a hybrid FV/FD scheme proposed by the authors [27,29] to simulate wave propagation through coastal porous structures. The present work is also regarded as a functional extension of our previous BT model developed for open-water problems [27], since the equations are roughly the same. A wide range of numerical simulations are carried out to test the capability of the Boussinesq equations, as well as the adopted numerical scheme. To this end, new experiments for solitary wave and regular wave interaction with a porous wall are also conducted in a wave flume.
The simulation of dam break flow in a porous abutment shows the merit of the proposed shock-capturing scheme. BT models based on the FD method do not resolve the initial discontinuity state, let alone the subsequent evolution process. A series of vertical two-dimensional and horizontal two-dimensional test cases for solitary wave and regular wave are further simulated. The computed results are found to be in generally good agreement with the experimental data, even compared with the RANS type simulation results in the literature. We thus conclude that the developed BT model is capable of handling wave interaction with porous structures in the coastal region for a wide range of scenarios after careful calibration of the porosity-related coefficients. GPU parallel computation further reduces the computation time greatly and enhances model capabilities in practical applications.