A Comparative Study of Numerical Models for Wave Propagation and Setup on Steep Coral Reefs

Complex factors including steep slopes, intense wave breaking, large bottom friction and remarkable wave setup should be considered while studying wave propagation over coral reefs, and how to simulate wave propagation and setup on coral reefs efficiently has become a primary focus. Several wave models can be used on coral reefs as have been published, but further testing and comparison of the reliability and applicability of these models are needed. A comparative study of four numerical wave models (i.e., FUNWAVE-TVD, Coulwave, NHWAVE and ZZL18) is carried out in this paper. These models’ governing equations and numerical methods are compared and analyzed firstly to obtain their differences and connections; then the simulation effects of the four wave models are tested in four representative laboratory experiments. The results show that all four models can reasonably predict the spectrum transformation. Coulwave, NHWAVE and ZZL18 can predict the wave height variation more accurately; Coulwave and FUNWAVE-TVD tend to underestimate wave setup on the reef top induced by spilling breaker, while NHWAVE and ZZL18 can predict wave setup relatively accurately for all types of breakers; NHWAVE and ZZL18 can predict wave reflection by steep reef slope more accurately. This study can provide evidence for choosing suitable models for practical engineering or establishing new models.


Introduction
The specific bathymetry makes the associated hydrodynamics on coral reefs more complex than that on normal coastal beaches (Yao et al., 2012). A typical fringing reef is characterized by a sloping reef face and a relatively gentle reef flat, as waves propagate from the deep sea to reef flats, water depth decreases drastically along reef slopes, inducing wave shoaling, reflection, refraction and nonlinear wave breaking. Wave breaking is accompanied by energy dissipation and can significantly change the wave energy spectrum, along with the flow changes. Most of the incident wave energy is dissipated through wave breaking and bottom friction; the breaking of waves also induces significant wave setup on reef flats. The wave setup is usually noticeable, especially at low tide, the wave setup can even be several times the original water depth on the reef top (Gourlay, 1996); conversely, the significant setup can change wave height. Both wave height and wave setup have a significant effect on the design elevation of the reef top structure.
With the demands for construction and protection of coral reefs, together with the trends of rising sea level, coral reef flattening and decreasing reef roughness, the study of the hydrodynamics over coral reefs has been a hot issue. Many scholars have conducted field observations (Zhu et al., 2004), laboratory experiments (Buckley et al., 2015;Demirbilek et al., 2007), and numerical simulations. Laboratory experiment can provide reliable results but with a highcost. It is also limited by the effect of scale and cannot be used under extreme conditions. Under this background, evaluating the reliability and applicability of existing numerical wave models can provide a reference for choosing suitable models to simulate wave propagation and setup on coral reefs, which is of great significance to practical engineering. Numerical wave models include phase-averaged models and phase-resolving models. Phase-averaged models simulate waves in a stochastic manner and have high efficiency, but they cannot predict infragravity waves and result in poor prediction of wave setup (Buckley et al., 2014). As a result, phase-resolving models are better choices. Navier-Stokes models are the most advanced phase-resolving models, and can be divided into two categories: grid-based methods and meshfree methods. Both categories can simulate splash-up and plunging due to wave breaking. Grid-based Navier-Stokes models need to employ free surface tracking techniques while in meshfree models plenty of nodes are required in the simulation domain to guarantee the simulation accuracy, causing the Navier-Stokes models too computational intensive to apply on field scale reef profiles. The published Navier-Stokes modeling studies on waves over coral reefs mainly target at the detailed transformation process of a solitary wave (Kazolea et al., 2016).
Among all the phase-resolving models, Boussinesq-type wave models have become the most prevailing wave models for coral reef systems due to their high efficiency and effectiveness. The development of Boussinesq equations suitable for steep slopes enhances their applicability on coral reef terrains (Zhu and Hong, 2001;Zou and Fang, 2008). Since the first application of Boussinesq equations on coral reefs by Skotner and Apelt (1999), many scholars have employed Boussinesq-type wave models to simulate waves over a coral reef system. Special attention should be paid to the governing equations, numerical schemes and wave breaking treatment methods because of the special terrains of coral reef systems. As for governing equations, Boussinesq equations with improved nonlinearity or Boussinesq equations more suitable for rapidly varying terrains were employed to build numerical models (Fang et al., 2016;Yao et al., 2012;Zhang et al., 2018); as for the numerical scheme, the shock-capturing finite volume (FV) method has been widely used recently, which improves the numerical stability and simulation precision (Roeber and Cheung, 2012); as for wave breaking treatment, a steep slope changes both the surf zone and the breaking threshold, and the breaking models were modified by adjusting the breaking parameters accordingly (Metallinos et al., 2016).
Non-hydrostatic models are proceeding rapidly and have been applied on simulating hydrodynamics of coral reef systems. Non-hydrostatic models assume the free surface to be a single value function of the horizontal coordinates, therefore the wave surface can be determined by applying the free surface boundary condition instead of using the free surface tracking techniques, and non-hydrostatic models are computationally more efficient than the Navier-Stokes models (Ma et al., 2012). Two representative open source non-hydrostatic wave models SWASH and NHWAVE have been applied in laboratory simulations and field simulations. Between the two models, SWASH requires a comparatively high vertical resolution (10-20 layers) to accurately simulate wave breaking, otherwise, the initiation of wave breaking is often delayed and wave energy dissipation in the surf zone is underestimated; while NHWAVE allows the model to accurately simulate wave breaking using fewer vertical layers owing to its more advanced numerical model scheme and grid system (Ma et al., 2012), which is an essential reason for choosing NHWAVE as a comparison model in this study.
Comparative studies of the applicability and prediction effect of various numerical wave models on waves over coral reefs are still limited. Kazolea et al. (2016) compared five numerical wave models by simulating a solitary wave breaking and overtopping over a fringing reef; the simulation results of all the models show a reasonable agreement with the measured data and there is no greater advantage from the Navier-Stokes models compared with the Boussinesqtype models. Buckley et al. (2014) evaluated the applicability of three numerical wave models (SWASH, SWAN, and XBeach) to coral reefs, the simulation results showed that the phase-averaged models can simulate the wave height variation reasonably after tuning the wave-breaking parameters; however, they cannot predict spectral transformation and infragravity wave motions, leading to poor simulation effect of the wave setup. Sheremet et al. (2011) compared the phase-averaged and phase-resolving mild-slope models and pointed out that the phase-resolving model predicted better results.
Previous comparative studies on wave models are often based on a single test case, causing them to unfavorably find factors influencing the prediction effect. In addition, the comparisons are often based on a single example, and thus cannot evaluate the applicability of numerical models on different bathymetries or different breaker types. In this paper, four wave models are selected according to the characteristics of different models, including three Boussinesqtype models and a non-hydrostatic model. The four models selected are FUNWAVE-TVD, Coulwave, NHWAVE and the numerical model established by Zhang et al. (2018) (referred to as ZZL18 hereafter). The governing equations, the numerical schemes and wave breaking treatment methods of selected models are compared separately, and the differences and connections between the models are summarized accordingly; and then the four models are applied to four typical test cases and their reliabilities are analyzed by comparing the simulated results and laboratory measurements. Combining the theoretical analysis of the numerical models, we can further obtain the factors influencing simulation effects. The conclusion of this study provides guiding significance for selecting appropriate model in practical engineering or developing new models suitable for simulating hydrodynamics on coral reefs.

Comparative analysis of the wave models' governing equations and numerical methods
A comprehensive comparative analysis of the four models' governing equations and numerical methods is presented in this section. Numerical methods include the numerical schemes, wave breaking treatment methods and boundary conditions. The description of the four selected models is omitted here for space limitation, readers can refer to Shi et al. (2012a), Ma et al. (2014), Lynett et al. (2008) and Zhang et al. (2018) for further details. After detailed subentry comparison of the models in Sections 2.1-2.4, the main differences and connections between them are summarized in Section 2.5.

Governing equations
Among the four selected numerical wave models, FUN-WAVE-TVD, Coulwave and ZZL18 use Boussinesq equations as the governing equations (one-layer mode and multilayer mode are included in Coulwave, and one-layer mode is used in this paper), and Boussinesq-type models are depth-intergraded two-dimensional models. The one-dimensional cases have included the primary wave transformation process (i.e., shoaling, breaking, reflection, generation of higher harmonics and infragravity waves), and only one-dimensional cases are considered for simplicity and clarity in this paper.
In one-dimensional cases, both the governing equations of FUNWAVE-TVD and Coulwave recover the equations of Wei et al. (1995), a moving reference level is employed in both models as in Kennedy et al. (2001); ZZL18 takes the Boussinesq equations of Kim et al. (2009) as the governing equations (referred to as KLS09 equations hereafter). The one-dimensional conservative form of Wei et al.'s (1995) equations reads η t + P x = 0; (1) is the horizontal mass flux; is the total water depth, is the still water depth and is the surface elevation; is the velocity at a reference , with , , ; is the depth-dependent second-order correction; and are both dispersive terms; for their detailed expressions, refer to Shi et al. (2012a).
KLS09 equations are suitable for rapidly varying bathymetry by including both the bottom curvature and squared bottom slope terms, the one-dimensional conservative form of the KLS09 equations reads where is the surface elevation, is the still water depth, represents the total water depth; is the horizontal mass flux; is the depth-averaged velocity and is the gravitational constant. Subscripts and are partial derivatives; is a free parameter of optimized dispersion with the value ; is the source term of dispersion. The comparison results of Wei et al.'s (1995) equations and KLS09 are as follows.
Both sets of equations achieve a Pade [2, 2] approximation to the exact dispersion relation, and are applicable to a water depth of with reasonable accuracy, where represents incident wavelength. The shoaling properties of the two sets of equations are both accurate to , the practical depth limitation of KLS09 is , while Wei et al.'s (1995) practical depth limitation is . Fig. 1 compares the wave speed and shoaling coefficient of the two sets of equations. Wei et al.'s (1995) equations retain terms of , and are second-order full-nonlinear equations; KLS09 equations retains terms of , and are weakly nonlinear equations. KLS09 equations are extended Boussinesq equations, in which both the bottom curvature and squared bottom slope terms are included. Numerical results show that KLS09 equations can predict wave reflection more accurately on very steep slopes (Kim et al., 2009).
NHWAVE takes the incompressible Navier-Stokes equations as the governing equations. Since the nonlinearity and dispersion properties of non-hydrostatic wave models still lack theoretical analysis at present, Ma et al. (2012)   proved by numerical experiment that NHWAVE can reflect the nonlinearity and dispersion properties of the wave with three vertical layers. With surface elevation errors restricted to 5%, the practical depth limitation of NHWAVE with three vertical layers can be .

Numerical Schemes
For spatial discretization, the Boussinesq-type models employ different schemes for the convection terms and the dispersion terms. The convection terms are discretized using a high-order shock-capturing FV method for all three Boussinesq-type models, a fourth-order MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) scheme combined with a HLL approximate Riemann solver is adopted to achieve the spatial scheme. A detailed description of the MUSCL scheme and HLL method can be found in Shi et al. (2012a).
The dispersion terms are discretized by a second-order center difference method in ZZL18 and FUNWAVE-TVD; while in Coulwave, the cell-averaged finite volume method is implemented and the second-order accuracy is selected, which is equivalent to the second-order center difference method. That is, the spatial discretization schemes of all the Boussinesq-type models' dispersion terms are consistent.
NHWAVE also adopts a hybrid FV/FD scheme for spatial discretization but of a lower-order. The second-order Godunov-type linear construction method along with the HLL Riemann approximation method is adopted to estimate the fluxes. The Poisson equation, which represents the non-hydrostatic pressure part, is discretized using a secondorder center finite-difference method.
The time integration schemes used include an implicit predictor-corrector scheme and explicit Runge-Kutta scheme; numerical tests verify that the third-order Runge-Kutta scheme can fully replace the predictor-corrector scheme (Shi et al., 2012b). Therefore, the time integration scheme of the three Boussinesq-type models can be considered to be consistent. Both the spatial discretization and time integration scheme of NHWAVE are of lower-order and finer grids are needed to achieve the same accuracy.

Wave breaking treatment
Wave breaking is a significant wave transformation process over coral reefs; in the strongest nonlinear surf zone, the nonlinearity lies in the energy dissipation from wave breaking. Energy dissipation induces wave height variation and the radiative stress gradient induced by wave breaking is a significant factor influencing the wave-induced setup. As a result, whether a model can accurately simulate the wave breaking process reflects the model's applicability to the coral reef terrains. There are two evaluation indexes for the models' breaking simulation effect: (1) whether the models can capture the breaking point accurately; (2) whether the model can predict the energy dissipation by wave breaking.
There are mainly two types of breaking treatment methods in the Boussinesq-type models. The first type simulates the energy dissipation caused by wave breaking by adding an extra dissipation term in the momentum equation, such as in the roller models and eddy viscosity models (Kennedy et al., 2000). The second type can be called a hybrid wave breaking model, wherein the governing equations are discretized using shock capturing schemes and wave breaking is judged by a user-defined threshold. In the swash zone, the Boussinesq equations are switched to nonlinear shallow water equations and the breaking wave was described as a bore or hydraulic jump. The advantage of the hybrid wave breaking model is its relative simplicity and effectiveness, but with fewer adjustable parameters the simulation precision may decrease for different terrains.
Researchers have compared the simulation effect of the eddy viscosity model and the hybrid wave breaking model (Kazolea and Ricchiuto, 2018). The comparison results show that the type of Boussinesq equations used has a significant influence on the hybrid wave breaking model; what is more the hybrid model is very sensitive to the grid size and strong oscillations may appear with fine grids. The eddy viscosity model demonstrates more robustness and stability when applied to various forms of equations or different meshing selections, and provides more accuracy in the simulation of wave height and water flow after wave breaking.
Coulwave and ZZL18 employ the eddy viscosity model to simulate wave breaking; the one-dimensional form of eddy viscosity model's expression is given by where is the eddy viscosity and ; is the empirical coefficient governing the magnitude of wave breaking w; parameter controls the occurrence of energy dissipation with a smooth transition from 0 to 1. The expression for is given by where determines the onset and stoppage of the breaking process and is evaluated by where is the moment when wave breaking begins; is the duration of wave breaking; is the transition time; and represent the threshold values of the occurrence and cease of wave breaking, respectively. The breaking parameters should be calibrated for excellent simulation effect, especially when is a key parameter influencing the simulation accuracy. Another point to be noted is that, on the steep slopes, the parameter controlling the occurrence of breaking is very different from that on a gentle slope and can be modified according to the breaking type.
FUNWAVE-TVD employs the hybrid wave breaking model and uses surface elevation to wave depth ratio as the breaking criterion. The cells are labeled as breaking cells when the ratio is larger than a certain threshold and dispersive terms are deactivated. The default value is in FUNWAVE-TVD and the value decreases on steep slopes. The breaking parameters are determined by the default values or trial calculation, and the parameters for different cases in this paper are listed in Table 1. For random wave cases, the default parameters or trimming ones have met the precision requirement; while a significant adjustment to the parameters is required for regular cases.
The non-hydrostatic model NHWAVE employs the shock-capturing scheme for spatial discretization, allowing it to capture wave breaking without relying on empirical formulations. The non-hydrostatic model with a Godunov-type scheme can predict wave height distribution, turbulence and undertow under breaking waves at least as accurate as the VOF model (Ma et al., 2012).

Boundary conditions
The boundary conditions involved include the bottom friction, wave generation, wave absorption and wet-dry interface.
The surfaces of coral reefs usually have large bottom roughness and the bottom friction changes with spatial distribution. The bottom friction has a significant influence on the wave height distribution and wave setup, and its unified form reads is the coefficient of bottom friction and is the velocity. In Boussinesq-type models, is the depth-integrated velocity; while in NHWAVE, is the bottom layer velocity. For spatially varied bottom friction, a spatially uniform friction coefficient can be used for each model.
An internal source wave maker guarantees the mass conservation of water, which ensures the accurate simulation of the wave-induced setup on the reef flat; in addition, an internal wave maker with a sponge layer can eliminate the influence of wave reflection from the reef-face, so that regular and random waves are made with the internal source function method in this paper. All four selected models can treat the wet-dry interface, using the linear extrapolation method or the wetting-drying scheme.

Main differences and connections between different
wave models Through the detailed subentry comparison of the models, it can be concluded that NHWAVE are different from the three Boussinesq-type models in several aspects; while the three Boussinesq-type models can be considered to have a single factor difference between each pair: the difference between FUNWAVE-TVD and Coulwave lies in their wave breaking treatment; the difference between ZZL18 and Coulwave lies in the governing equations. The main characteristics of different wave models are listed in Table 2 for clarity.

Comparison of numerical models on wave propagation over fringing reefs
Four typical laboratory experiments are utilized to comprehensively test and compare the performance of different models on coral reefs. When the fore-reef slope is steeper than 1:1, the coral reef is considered very steep. According to the characteristics of incident waves and topography, the four test cases are termed regular wave propagation over steep fringing reefs, regular wave propagation over very steep fringing reefs, random wave propagation over composite-slope fringing reefs and random wave propagation on steep fringing reefs.

Regular wave propagation over steep fringing reefs
A representative case for regular wave transformation on an idealized fringing reef by Yao et al. (2012) is simulated. The reef profile and positions of the wave gauges are shown in Fig. 2. Twelve wave gauges were positioned and marked WG1 to WG12 from left to right. The reef slope was 1:6, the off-reef water depth m, the water depth over the reef flat was 0.1 m, and sponge layers were installed at both sides of the domain. The incident wave height was 0.095 m and the wave period was 1.25 s. dx = 0.03 A grid size of m is selected in all the models, and three vertical layers are used in NHWAVE. The models were run for 200 wave periods and the last 125 wave periods were used for data analysis. The computed wave heights and mean water level (MWL) of all the models are compared with the experimental data in Fig. 3a and Fig. 3b, respectively.
As shown in Fig. 3, all the models' simulated wave heights and locations of breaking point agree well with laboratory data; significant standing waves appear in front of the slope due to the reflection; in the surf zone, each model's simulation precision decreases: the observed wave height decreases rapidly after breaking, while the simulated results underestimate the speed of decrease in the wave height; all the models underestimate the stable wave height on the reef flat and ZZL18 has the highest precision.
Regarding MWL, all the models' simulation results are almost identical in the off-reef deep water; but obvious differences appear in the surf zone, where the MWL simulation precision of the NHWAVE model and the Coulwave model are higher. All models underestimate the maximum set-down at the breaking point; on the reef flat, FUNWAVE-TVD underestimates the wave setup, and the other models' simulation results are similar, among them ZZL18 predicts the maximum wave setup with the highest accuracy.
To quantify the performance of different models, the model skill parameter is introduced (Gallagher et al., 1998). For a nearly perfect model result, the value of model skill should approach 1; the expression of model skill is X meas X pred N where represents the measured value, represents the calculated value, and is the total number of observations. There were twelve wave gauges in Yao et al.'s (2012) experiment, labeled as WG1 to WG12 from the left to right with locations as shown in Fig. 2. WG5 to WG8 are located in the surf zone, while WG1 to WG4 are located before wave breaking point and WG9 to WG12 are located behind the surf zone (Yao et al., 2012). Comparing the wave height and wave setup in three zones separately is helpful to clarify the performance of the models. Table 3 and Table 4   ZHANG Shan-ju et al. China Ocean Eng., 2019, Vol. 33, No. 4, P. 424-435 pare the model skills of wave heights and wave setup in three different zones. Data analysis obtains the same conclusion as the observation of the figure: the wave height prediction precisions of all the models are relatively high before and after wave breaking with model skill values above 0.8; in the surf zone, the precisions are lower since the models underestimate the wave height decrease speed. The wave setup prediction precision of FUNWAVE-TVD is significantly lower than that of the other three models. The wave setup predicted by NHWAVE and that of ZZL18 agrees best with observations on the reef.

Regular wave propagation over very steep fringing reefs
Wen et al. (2019) conducted laboratory experiments of regular wave on idealized coral reefs, and a significant case is modelled here. Fig. 4 shows the diagram of the experimental terrain and the positions of the wave gauges. The off-reef water depth was 1.9 m, and the reef height was 1.8 m and was followed by a lagoon. The reef slope was 1:1 and the reef-flat is 8.6 m long. The incident waves were regular waves with wave height and wave period .
In this laboratory experiment, the lagoon is connected with the deep sea before the reef slope, so the wave setup would not induce significant wave set-down before the reef slope, and one dimensional simulation cannot deal with this problem well. To solve this problem, the periodic boundary conditions were employed by Wen et al. (2019). This problem is proceeded using a simple treatment in this paper: the initial water level was finely tuned to guarantee the water level off-reef to be consistent with the measurement after making waves. The leeside vertical faces of the reef were modified to a 1:1 slope to ensure numerical stabilities. A grid size of m is selected in the Boussinesq-type models; while a grid size of m and three vertical layers are used in NHWAVE. The wave heights and wave setup are presented in Fig. 5a and Fig. 5b. Fig. 5a shows that all the models can predict the location of wave breaking accurately, and the wave height variation is also captured reasonably. NHWAVE and ZZL18 can predict wave reflection before the reef more accurately than the other two models; the breaking wave height predicted by NHWAVE and FUNWAVE-TVD are bigger; on the reef flat, FUNWAVE-TVD underestimate the wave height. Fig. 5b shows that all the models' simulated MWL are consistent with measurements overall; the maximum setdown by NHWAVE and FUNWAVE-TVD are more obvious; NHWAVE and ZZL18 overestimate the wave setup on reef flat. Table 5 compares the model skill for predicting the wave height and wave setup for different models. All the models' predicted wave heights are relatively accurate. NHWAVE and ZZL18 are more accurate because they can better predict wave reflection before the reef and the wave height after breaking. All the models can predict wave setup with reasonable precision, of which Coulwave performs better, as it better predicted the wave setup on reef flat.

Random wave propagation over composite-slope fringing reefs
0.031 m Demirbilek et al. (2007) reported a laboratory experiment on transformation of random waves over a fringing reef, and the corresponding bottom topography is shown in Fig. 6. Nine wave gauges were positioned and marked as WG1 to WG12. The right boundary was a 1:12 slope. The off-reef water depth was 0.531 m, and the topography transitioned from a composite reef face with an average slope of 1:10.2 to a flat reef top with a water depth of . The significant wave height of incident random waves was 0.075 m and the peak wave period was 1.5 s.
Fourier analysis is conducted on the time series of the wave surface at WG1 to generate the corresponding intern-   al wave signal for the desired wave. A grid size of m is selected in NHWAVE, while the other three Boussinesqtype models employ a grid size of . The present simulation lasted for and the results from the last were used for data analysis. The computed and measured significant wave heights and MWL are shown in Fig. 7a and Fig. 7b, respectively. The measured surface elevations at WG4 appeared to be faulty and are therefore not used here (Zijlema, 2012).
From the preceding figures, we can see that all models' computed wave heights and wave setup agree well with the measurements. Especially for the wave heights, all the models' simulated wave heights are almost identical in the surf zone and at the reef flat, and only a tiny distinction exists in the fore-reef zone. Every model can predict the wave setup on the reef flat accurately; Coulwave fails to predict the setdown before wave-breaking; the other three models' simulation results embody the overall trend of set-down before wave-breaking and setup on the reef flat. Table 6 shows the model skill values for predicting the significant wave heights and wave setup, the wave setup on the reef flat is listed out separately because it is the main concern for practical engineering. The simulated wave heights from all the models are highly comparable with measurements. All the model skill values are around 0.95, and Coulwave's model skill is the highest at 0.967. The four models' simulation precisions of wave setup are lower than that for wave heights; the over model skill values of ZZL18 and NHWAVE are the highest, and all the models' prediction precisions on the reef flat are rather high.
The breaking of random waves over coral reefs is usually accompanied by energy transfer. Fig. 8 compares the measured and simulated wave energy spectra of four representative wave gauges (WG3, WG6, WG8, and WG9). As waves propagate from deep water to reef flats, wave energy gradually transfers from the peak frequency to both higher and lower frequencies. Waves break near the reef edge, and most of the wave energy around the peak frequency of incident waves is dissipated in the surf zone area. The wave energy is dominated by infragravity waves on the reef flat. All the models can well simulate the wave energy transfer in this process.
To test the models' simulation effects on wave transformation, the surface elevations of two representative wave gauges are compared. Fig. 9 and Fig. 10 compare the time series of all the models' computed surface elevations from 100 s to 150 s with the measured data at WG6 and WG8. η > 2 cm Fig. 9 shows that all the models can simulate the trend towards bigger wave height and greater wave steepness due to nonlinear shoaling. Both the predicted wave shapes and the wave phases agree well with the measurements. Comparing the predicted surface elevations with the observations in Fig. 9, we can see that when the surface elevation is low, the predicted results agree well with the observations; the discrepancy begins to increase when the surface elevation increases to , probably because some of the waves with high surface elevations start to break. WG8 is located on the reef flat after wave-breaking; in this area, the orders of the predicted values and observations are close, but relatively large discrepancies exist in the wave shapes and wave phases, as shown in Fig. 10.

Random wave propagation over steep fringing reefs
A typical case in Smith et al.'s (2012) experiment is selected for the simulation, corresponding to a 0.439 m offshore water depth, a 1:5 reef slope and a rough surface. The reef flat is 7.3 m long with varying gentle slopes, and the boundary at the far right is a 1:10 slope. Fig. 11 shows a diagram of the experimental terrain and the positions of the wave gauges. A total of 12 wave gauges were positioned and labeled WG1 to WG12 from left to right. Random waves we-   ZHANG Shan-ju et al. China Ocean Eng., 2019, Vol. 33, No. 4, P. 424-435 re generated based on a TMA spectrum with a significant wave height of 0.0749 m and a peak wave period of 1.41 s. Fourier analysis is conducted on the time series of the wave surface at WG1 to generate the corresponding intern- al wave signal for the desired wave. A grid size of and a time step of s are selected in all the models, and different bottom friction coefficients are employed in the off-reef ( ) and reef flat ( ) areas. The computed and measured significant wave heights and MWL are shown in Fig. 12a and Fig. 12b, respectively. The model skill values for the predicted significant wave height and wave setup are shown in Table 7.
From the preceding figures, we can see that the predicted significant wave heights from different models all Table 7 Model skill comparison of significant wave height and wave setup for different models in Smith et al.'s (2012)       agree well with the experimental data. The model skill values of all the models' predicted significant wave heights are larger than 0.9, and the prediction precision of NHWAVE is slightly higher than that of the other models on the reef flat.
As for the wave setup on the reef flat, ZZL18 and NHWAVE can simulate the wave setup with relatively high accuracy, while Coulwave and FUNWAVE-TVD underestimate the wave setup on the reef flat. The wave setup has a certain effect on the wave heights on the reef flats. Especially for the reformed stale waves at WG12, there is a positive correlation between the wave setup and wave height. Every model can predict the energy transfer and the generation of infragravity waves in this case; the results can be seen in Su et al. (2015), Zhang et al. (2018) and Ma et al. (2014) and therefore they are omitted here.

Discussion
h/L 0 < 0.5 The applicability and prediction precision of four numerical wave models are evaluated only by laboratory experiments in this study. When applied on the actual filed terrains, the nonlinear-dispersive properties of Coulwave and NHWAVE can be improved for deep water by choosing multilayer mode or by increasing vertical layers. The dispersion properties of ZZL18 and FUNWAVE-TVD are applicable to a water depth of , when off-reef water is deeper than half-wavelength, topography variation has little effect on wave speed, and the water depth can be taken as half-wavelength. By doing this, the dispersion properties of ZZL18 and FUNWAVE-TVD meets accuracy requirements for regular waves in deep water, as to random wave, the simulation precision decreases due to the errors of high frequency waves. The nonlinear properties of both FUNWAVE-TVD and ZZL18 are not accurate in deep water. This study concerns on the wave transformation on steep fore-reef, the disadvantage that FUNWAVE-TVD and ZZL18 lack of nonlinear-dispersive precisions in deep waters can be solved by nesting with other models (such as multilayer Coulwave), both FUNWAVE-TVD and ZZL18 provide interfaces for nested models.
After solving the depth limitation problem of the wave models, the applicability of ZZL18 on fore-reef slopes should be further discussed because its governing equations are weakly nonlinear. Significant differences exist between the characteristics of steep slope and gentle slope wave-breaking, the horizontal distance for wave transformation is very short and wave breaking begins before the wave height increases adequately. For example, the horizontal distance for wave transformation before wave breaking in Wen et al.'s (2019) experiment is only 0.3 times the incident wavelength. As a result, the strong nonlinearity of the wave lies in depth-induced wave breaking on steep slopes. Wave breaking is treated using energy dissipation methods in Boussinesq models and cannot simulate the splash-up and plunging in the breaking process. Therefore, the nonlinearity of the governing equation has little influence, while simulating wave reflection accurately and treating wave breaking reasonably are more important. Previous researches have pointed out that the inclusion of higher nonlinear terms in Boussinesq equations does not necessarily improve the MWL (Yao et al., 2012). The reflection increases as the slope increases; for the slopes steeper than 1:1, ZZL18 can simulate wave reflection more accurately (Kim et al., 2009), therefore it has higher simulation precision to non-breaking cases. Therefore, ZZL18 is applicable to coral reefs though its governing equations are weakly nonlinear.
When the basic type of wave breaking is surging breaker or spilling breaker, all the models can predict wave setup with reasonable accuracy. When the basic wave breaking type is plunging breaker, significant differences exist between different models: FUNWAVE-TVD tends to underestimate the wave setup after wave breaking. This problem was also found by other researchers and was widely believed that this could be a deficiency of the hybrid wave breaking model; the underestimate of wave setup is also found in Coulwave; ZZL18 and NHWAVE can predict wave setup with relatively high accuracy for different breaker types. Limited by its weakly nonlinearity, ZZL18 cannot predict the relatively low-frequency waves as accurately as other models. Based on the comparison of predicting effect of different models as well as their stability, we can conclude that NHWAVE and ZZL18 behave better than Coulwave and FUNWAVE-TVD, and Coulwave behaves better than FUN-WAVE-TVD. Among the three Boussinesq-type models, the main difference between ZZL18 and Coulwave lies in their governing equations, in which we conclude that the governing equations of ZZL18 are more suitable for steep bathymetry; the main difference between Coulwave and FUNWAVE-TVD lies in their wave breaking treatment methods, in which we conclude that an eddy viscosity model simulates the wave setup more accurately than a hybrid wave breaking model on steep slopes.
On steep slopes, the breaking of random wave can still be reasonably simulated using the default breaking parameter; while the breaking of regular wave cannot be reasonably captured unless the breaking parameters are well adjusted. The reason for this difference is that regular waves break at a certain point, while the surf zone of random waves is an area, so regular waves require higher precision to capture wave-breaking compared with random waves. The initiation of wave breaking is related to the slope of the topography, and waves are easier to break on steeper slopes, and smaller value of or should be taken in eddy viscosity model or hybrid wave breaking model. NHWAVE can simulate wave breaking and wave transformation with reasonable precision by using three vertical layers, and the simulating effect of NHWAVE is stable under different wave conditions and different topographies. The three-dimensional NHWAVE has no obvious advant-age over the two-dimensional Boussinesq-type models on predicting wave height variation, wave setup and the spectrum transformation, but has a higher vertical resolution. NHWAVE needs finer grids to achieve the same precision as other Boussinesq-type models in the horizontal direction; what is more, solving the Poisson equation is the most timeconsuming in NHWAVE, making the computational cost more expensive than Boussinesq-type models.

Conclusions
In this paper, the applicability and prediction precision of four numerical wave models in complex coral reef environments are tested and compared by using four representative laboratory experiments. Wave reflection, wave shoaling, nonlinear interactions, energy dissipation by wave breaking and wave-induced setup are considered in the tests. By comparing the simulation results and the laboratory data, the following main conclusions can be drawn.
(1) Waves are easier to break on steep slopes; by selecting appropriate breaking parameter, all the selected wave models can predict the wave height variation and the spectrum transformation with reasonable accuracy.
(2) ZZL18 and NHWAVE can predict wave setup reasonably for different breaker types; Coulwave and FUN-WAVE-TVD can predict wave setup reasonably by spilling breaker and surging breaker, but tend to underestimate wave setup by plunging breaker. The predict precision of relatively low-frequency wave by ZZL18 is slightly lower than that of other models.
(3) If only the horizontally-varying spatial information such as variation of wave height and wave setup are considered, NHWAVE has no obvious advantage over the Boussinesq-type models, but the overall computation cost of NHWAVE is about one order of magnitude larger. However, NHWAVE has a higher vertical resolution and is therefore more suitable for cases considering the vertical flow structure.
(4) To the Boussinesq-type models, the breaking parameters should be carefully adjusted for regular waves breaking on steep slopes; while the breaking of random waves can be simulated reasonably using the default breaking parameters. Boussinesq wave models with eddy viscosity model can simulate wave propagation on steep coral reefs more accurately than those using hybrid wave breaking models.