Large eddy simulation modeling of tsunami-like solitary wave processes over fringing reefs

Many low-lying tropical and subtropical reef-fringed coasts are vulnerable to inundation during tsunami events. Hence accurate prediction of tsunami wave transformation and run-up over such reefs is a primary concern in the coastal management of hazard mitigation. To overcome the deficiencies of using depth-integrated models in modeling tsunami-like solitary waves interacting with fringing reefs, a three-dimensional (3-D) numerical wave tank based on the computational fluid dynamics (CFD) tool OpenFOAM^® is developed in this study. The Navier–Stokes equations for two-phase incompressible flow are solved, using the large eddy simulation (LES) method for turbulence closure and the volume-of-fluid (VOF) method for tracking the free surface. The adopted model is firstly validated by two existing laboratory experiments with various wave conditions and reef configurations. The model is then applied to examine the impacts of varying reef morphologies (fore-reef slope, back-reef slope, lagoon width, reef-crest width) on the solitary wave run-up. The current and vortex evolutions associated with the breaking solitary wave around both the reef crest and the lagoon are also addressed via the numerical simulations.

1 Introduction

A tsunami is an extremely destructive natural disaster, which can be generated by earthquakes, landslides, volcanic eruptions and meteorite impacts. Tsunami damage occurs mostly in the coastal areas where tsunami waves run up or run down the beach, overtop or ruin the coastal structures, and inundate the coastal towns and villages (Yao et al., 2015). Some tropic and subtropical coastal areas vulnerable to tsunami hazards are surrounded by coral reefs, especially those in the Pacific and Indian oceans. Among various coral reefs, fringing reefs are the most common type. A typical cross-shore fringing reef profile can be characterized by a steep offshore fore-reef slope and an inshore shallow reef flat (Gourlay, 1996). There is also possibly a reef crest lying at the reef edge (e.g., Hench et al., 2008) and/or a narrow shallow lagoon existing behind the reef flat (e.g., Lowe et al., 2009a). Over decades, fringing reefs have been well known to be able to shelter low-lying coastal areas from flood hazards associated with storms and high surf events (e.g., Cheriton et al. 2016; Lowe et al., 2005; Lugo-Fernandez et al., 1998; Péquignet et al., 2011; Young, 1989). However, after the 2004 Indian Ocean tsunami, the positive role of coral reefs in mitigating the tsunami waves has attracted the attention of scholars who conduct postdisaster surveys (e.g., Chatenoux and Peduzzi, 2007; Ford et al., 2014; Mcadoo et al., 2011). There is consensus among the researchers that, in addition to the establishment of the global tsunami warning system, the cultivation of coastal vegetation (mangrove forest, coral reef, etc.) is also one of the coastal defensive measures against the tsunami waves (e.g., Dahdouh-Guebas et al., 2006; Danielsen et al., 2005; Mcadoo et al., 2011). Numerical models have been proven to be powerful tools to investigate tsunami wave interaction with the mangrove forests (e.g., Huang et al., 2011; Maza et al., 2015; Tang et al., 2013, and many others). Comparatively speaking, their applications in modeling coral reefs subjected to tsunami waves are still very few.

Over decades, modeling wave processes over reef profiles face several challenges such as steep fore-reef slope, complex reef morphology and spatially varied surface roughness. Local but strong turbulence due to wave breaking in the vicinity of the reef edge needs to be resolved. Among various approaches for modeling wave dynamics over reefs, two groups of models are the most pervasive. The first group focuses on using the phase-averaged wave models and the nonlinear shallow water equations to model the waves and the flows, respectively, in field reef environments, and typically the concept of radiation stress (Longuet-Higgins and Stewart, 1964) or vortex force (Craik and Leibovich, 1976) is used to couple the waves and the flows (e.g., Douillet et al., 2001; Kraines et al., 1998; Lowe et al., 2009b, 2010; Van Dongeren et al., 2013; Quataert et al., 2015). As for modeling tsunami waves at a field scale, we are only aware of Kunkel et al.'s (2006) implementation of a nonlinear shallow water model to study the effects of wave forcing and reef morphology variations on the wave run-up. However, their numerical model was not verified by any field observations. The second group aims at using the computationally efficient and phase-resolving model based on the Boussinesq equations. This depth-integrated modeling approach employs a polynomial approximation to the vertical profile of velocity field, thereby reducing the dimensions of a three-dimensional problem by one. It is able to account for both nonlinear and dispersive effects at intermediate water level. At a laboratory scale, Boussinesq models combined with some semiempirical breaking wave and bottom friction models have been proven to be able to simulate the motions of regular waves (Skotner and Apelt, 1999; Yao et al., 2012), irregular waves (Nwogu and Demirbilek, 2010; Yao et al., 2016, 2019) and infragravity waves (Su et al., 2015; Su and Ma, 2018) over fringing reef profiles.

The solitary wave has been employed in many laboratory/numerical studies to model the leading wave of a tsunami. Compared to the aforementioned regular/irregular waves, the numerical investigations of solitary wave interaction with the laboratory reef profile are much fewer. Roeber and Cheung (2012) was the pioneer study to simulate the solitary wave transformation over a fringing reef using a Boussinesq model. Laboratory measurements of the cross-shore wave height and current across the reef as conducted by Roeber (2010) were reproduced by their model. More recently, Yao et al. (2018) also validated a Boussinesq model based on their laboratory experiments to assess the impacts of reef morphologic variations (fore-reef slope, back-reef slope, reef-flat width, reef-crest width) on the solitary wave run-up over the back-reef beach. Despite the above applications, several disadvantages still exist in using the Boussinesq-type models: (1) Boussinesq equations are subjected to the mild-slope assumption – thus using them for reefs with a steep fore-reef slope is questionable, particularly when there is a sharp reef crest located at the reef edge; (2) wave breaking could not be inherently captured by Boussinesq-type models – thus empirical breaking model or special numerical treatment is usually needed; (3) Boussinesq models could not resolve the vertical flow structure associated with the breaking waves due to the polynomial approximation to the vertical velocity profile.

To remedy the above deficiencies of using Boussinesq-type models to simulate the solitary processes (wave breaking, bore propagation and run-up) over the fringing reefs, we develop a 3-D numerical wave tank based on the computational fluid dynamics (CFD) tool OpenFOAM^® (Open Field Operation and Manipulation) in this study. OpenFOAM^® is a widely used open-source CFD code in the modern industry supporting two-phase incompressible flow (via its solver interFoam). With appropriate treatment of wave generation and absorption, it has been proven to be a powerful and efficient tool for exploring complicated nearshore wave dynamics (e.g., Higuera et al., 2013b). In this study, the Navier–Stokes equations for an incompressible fluid are solved. For the turbulence closure model, although large eddy simulation (LES) demands more computational resources than Reynolds-averaged Navier–Stokes (RANS) equations, it computes the large-scale unsteady motions explicitly. Importantly, it could provide more statistical information for the turbulence flows in which large-scale unsteadiness is significant (Pope, 2000). Thus the LES model is adopted by considering that the breaking-wave-driven flow around the reef edge/crest is fast and highly unsteady. The free-surface motions are tracked by the widely used volume-of-fluid (VOF) method.

In this study, we first validate the adopted model by the laboratory experiments of Roeber (2010) as well as our previous experiments (Yao et al., 2018). The robustness of the present model in reproducing such solitary wave processes as wave breaking near the reef edge/crest, turbulence bore propagating on the reef flat and wave run-up on the back-reef beach is demonstrated. The model is then applied to investigate the impacts of varying reef morphologies (fore-reef slope, back-reef slope, lagoon width, reef-crest width) on the solitary wave run-up. The flow and vorticity fields associated with the breaking solitary wave around the reef crest and the lagoon are also analyzed by the model results. The rest of this paper is organized as follows. The numerical model is firstly described in Sect. 2. It is then validated by the laboratory data from the literature as well as our data in Sect. 3. What follows in Sect. 4 are the model applications for which laboratory data are unavailable. The main conclusions drawn from this study are given in Sect. 5.

2 Numerical methods

2.1 Governing equations

To simulate breaking-wave processes across the reef, the LES approach is employed to balance the need of resolving a large portion of the turbulent flow energy in the domain while parameterizing the unresolved field with a subgrid closure in order to maintain a reasonable computational cost. The filtered Navier–Stokes equations are essential to separate the velocity field that contains the large-scale components, which is performed by filtering the velocity field (Leonard, 1975). The filtered velocity in the ith spatial coordinate is defined as

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{u}}_i\left(x\right)=\int G(x,x^{\prime })u_i\left(x^{\prime }\right)\mathrm{d}{x}^{\prime },\end{array}$$

where G(x, x^′) is the filter kernel, which is a localized function. The eddy sizes are identified using a characteristic length scale, Δ, which is defined as

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }=(\mathrm{\Delta }x\cdot \mathrm{\Delta }y\cdot \mathrm{\Delta }z{)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where Δx, Δy and Δz are the grid size in streamlines, spanwise and vertical directions, respectively. Eddies that are larger than Δ are roughly considered as large eddies, and they are directly solved. Those which are smaller than Δ are small eddies.

For incompressible flow, the filtered continuity and momentum equations are as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial x_i}}=\mathrm{0},\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\rho }{\stackrel{\mathrm{‾}}{u}}_i}{\partial t}}+{\displaystyle \frac{\partial \left(\mathit{\rho }{\stackrel{\mathrm{‾}}{u}}_i{\stackrel{\mathrm{‾}}{u}}_j\right)}{\partial x_i}}=-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial x_i}}+\mathit{\rho }{g}_i+\mathrm{2}\mathit{\mu }{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{S}}_{ij}}{\partial x_j}}-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}^r}{\partial x_j}},\end{array}$$

where ρ is the water density, μ is the dynamic viscosity, $\stackrel{\mathrm{‾}}{p}$ is the filtered pressure, ${\stackrel{\mathrm{‾}}{S}}_{ij}$ is the strain rate of the large scales defined as

$$\begin{array}{}\text{(5)}& {\stackrel{\mathrm{‾}}{S}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_j}{\partial x_i}}\right)\end{array}$$

and ${\mathit{\tau }}_{ij}^{\mathrm{r}}$ is the residual stress approximated by using subgrid-scale (SGS) models to get a full solution for the Navier–Stokes equations.

The residual stress is usually calculated by a linear relationship with the rate-of-strain tensor based on the Boussinesq hypothesis. The one-equation eddy viscosity mode, which is supposed to be better than the well-known Smagorinsky model for solving the highly complex flow and shear flow (Menon et al., 1996), is employed in the present study. Based on the one-equation model (Yoshizawa and Horiuti, 1985), the subgrid stresses are defined as

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{ij}^{\mathrm{r}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{k}_{\mathrm{S}}{\mathit{\delta }}_{ij}-\mathrm{2}{v}_t\left({\stackrel{\mathrm{‾}}{S}}_{ij}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\stackrel{\mathrm{‾}}{S}}_{kk}{\mathit{\delta }}_{ij}\right),\end{array}$$

where δ_ij is the Kronecker delta, and ν_t is the SGS eddy viscosity, which is given by

$$\begin{array}{}\text{(7)}& {\mathit{\nu }}_t=C_k\mathrm{\Delta }\sqrt{k_{\mathrm{S}}},\end{array}$$

and the SGS kinetic energy k_S needs to be solved by

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial k_{\mathrm{s}}}{\partial t}}+{\stackrel{\mathrm{‾}}{u}}_i{\displaystyle \frac{\partial k_{\mathrm{s}}}{\partial x_i}}={\displaystyle \frac{\partial }{\partial x_i}}\left({\displaystyle \frac{\mathit{\mu }}{P_{\mathrm{r}}}}{\displaystyle \frac{\partial k_{\mathrm{s}}}{\partial x_i}}\right)-{\displaystyle \frac{{\mathit{\tau }}_{ij}^{\mathrm{r}}}{\mathit{\rho }}}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_j}{\partial x_i}}-{\displaystyle \frac{C_{\mathit{\epsilon }}{k}_{\mathrm{s}}^{\mathrm{3}/\mathrm{2}}}{\mathrm{\Delta }}},\end{array}$$

where C_k=0.094, C_ε=0.916 and P_r=0.9 as suggested by the OpenFOAM^® user guide (OpenFOAM Foundation, 2013).

The presence of the free-surface interface between the air and water is treated through the commonly used VOF method (Hirt and Nichols, 1981), which introduces a volume fraction and solves an additional modeled transport equation for this quantity. The general representation of fluid density ρ is written as

$$\begin{array}{}\text{(9)}& \mathit{\rho }=\mathit{\alpha }{\mathit{\rho }}_{\mathrm{1}}+(\mathrm{1}-\mathit{\alpha }){\mathit{\rho }}_{\mathrm{2}},\end{array}$$

where ρ₁=1000 kg m⁻³ is the density of water, ρ₂=1 kg m⁻³ is the density of air and α is the volume fraction of water contained in a grid cell. The distribution of α is modeled by the advection equation

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial \mathit{\alpha }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\alpha }{\stackrel{\mathrm{‾}}{u}}_i\right)+\mathrm{\nabla }\cdot \left[\mathit{\alpha }(\mathrm{1}-\mathit{\alpha })u_i^{\mathrm{r}}\right]=\mathrm{0}.\end{array}$$

The last term on the left side is an artificial compression term, avoiding the excessive numerical diffusion and the interface smearing; the new introduced $u_i^{\mathrm{r}}$ is a velocity field suitable to compress the interface.

In the present solver interFoam, the algorithm PIMPLE, which is a mixture of the PISO (Pressure Implicit with Splitting of Operators) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithms, is employed to solve the coupling of velocity and pressure fields. The MULES (multidimensional universal limiter for explicit solution) method is used to maintain boundedness of the volume fraction independent of the underlying numerical scheme, mesh structure, etc. The Euler scheme is utilized for the time derivatives, the Gauss linear scheme is used for the gradient term and the Gauss linear corrected scheme is selected for the Laplacian term. The detailed implementation can be founded in the OpenFOAM^® user guide (OpenFOAM Foundation, 2013).

2.2 Wave generation and absorption

Wave generation and absorption are essentials for a numerical wave tank, but they are not included in the official version of OpenFOAM^®. Therefore, supplementary modules were developed by the other users, e.g., waves2Foam (Jacobsen et al., 2012) and IH-FOAM (Higuera et al., 2013a). In this study, the IH-FOAM is selected because it employs an active wave absorbing boundary and does not require an additional relaxation zone as used by waves2Foam. Meanwhile, it supports many wave theories including the solitary wave theory. The free surface and velocity for a solitary wave generation in IH-FOAM are (Lee et al., 1982)

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }\mathit{\eta }=H\sec h^{\mathrm{2}}\left(\sqrt{{\displaystyle \frac{\mathrm{3}H}{\mathrm{4}{h}^{\mathrm{3}}}}}X\right),\text{(12)}& {\displaystyle }& {\displaystyle \frac{u}{\sqrt{gh}}}={\displaystyle \frac{\mathit{\eta }}{h}}\left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\displaystyle \frac{\mathit{\eta }}{h}}+{\displaystyle \frac{h}{\mathrm{3}}}{\displaystyle \frac{h}{\mathit{\eta }}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\displaystyle \frac{z^{\mathrm{2}}}{h^{\mathrm{2}}}}\right){\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\eta }}{\mathrm{d}{X}^{\mathrm{2}}}}\right],\text{(13)}& {\displaystyle }& {\displaystyle \frac{w}{\sqrt{gh}}}={\displaystyle \frac{-z}{h}}\left[\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{h}}\right){\displaystyle \frac{\mathrm{d}\mathit{\eta }}{\mathrm{d}X}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{h}^{\mathrm{2}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{z^{\mathrm{2}}}{h^{\mathrm{2}}}}\right){\displaystyle \frac{{\mathrm{d}}^{\mathrm{3}}\mathit{\eta }}{\mathrm{d}{X}^{\mathrm{3}}}}\right],\end{array}$$

where η is the free-surface elevation, H is the wave height, h is the water depth, $X=x-ct$, $c=\sqrt{g(h+H)}$ is the wave celerity, and u and w are the velocities in the streamwise and vertical directions, respectively.

3 Model validation

3.1 Experimental settings

The first set of laboratory experiments serving as validation is from Roeber (2010), who reported two series of experiments conducted at Oregon State University, USA, in separate wave flumes. In this study, we only reproduce their experiments in the large wave flume, which is 104 m long, 3.66 m wide and 4.57 m high. As illustrated in Fig. 1a, the two-dimensional (2-D) reef model, starting at 25.9 m from the wavemaker, was built by a plane fore-reef slope attached to a horizontal reef flat 2.36 m high followed by a back-reef vertical wall. Both the waves and flows across the reef profile were measured by 14 wave gauges (wg1–wg14) and 5 ADVs (acoustic Doppler velocimeters), respectively. Only two scenarios for the reef with and without a trapezoidal reef crest subjected to two incident waves are reported in this study (see also Table 1). The large wave flume experiments enable us to test our model's ability to handle relatively large-scale nonlinear dispersive waves together with wave breaking, bore propagation and associated wave-driven flows. For a more detailed experimental setup, see Roeber (2010).

Figure 1 Experiment settings for (a) Roeber (2010) and (b) Yao et al. (2018).

Table 1 Reef configuration and wave condition for the tested scenarios.

The second set of 2-D reef experiments for model validation comes from our previous work (Yao et al., 2018). These experiments were conducted in a small wave flume 40 m long, 0.5 m wide and 0.8 m high at the Changsha University of Science and Technology, P. R. China. As shown in Fig. 1b, a plane slope was built 27.3 m from the wavemaker and it was truncated by a horizontal reef flat 0.35 m high. A back-reef beach of 1:6 was attached to the end of the reef flat. The surface elevations were measured at eight cross-shore locations (G1–G8) and no flow measurement was performed. However, a CCD camera was installed to record the process of water uprush on the back-reef slope. Thus the model's robustness to capture the whole process of solitary wave transformation over the reef flat and run-up on the back-reef beach can be evaluated. In this study, we only simulate the tested idealized reef profile with and without a lagoon at the rear of reef flat subjected to the same wave condition (see also Table 1). The lagoon was formed by two 1:1 slopes connecting the reef flat and the toe of the back-reef beach to the flume bottom. The dimensions of the fore-reef slope and the reef flat, the water depths over the reef flat and the incoming wave heights were designed according to the Froude similarity with a target geometric-scale factor of 1:20. See Yao et al. (2018) for the detailed laboratory settings.

3.2 Numerical settings

By considering a balance between the computational accuracy and efficiency, the computational domain (Fig. 2a) is designed to reproduce the main aspects of the laboratory settings. We calibrate the model with the principle that the computed leading solitary wave height at the most offshore gauge should exactly reproduce its measurement. For a solitary wave, wave length (L) can be estimated as a distance containing 95 % of the total mass of the solitary wave, which yields $L=\mathrm{2.12}h/\sqrt{H_i/h}$. The largest offshore wave length according to the wave conditions in Table 1 is L=8.44 and 1.52 m for the scenario of Roeber (2010) and Yao et al. (2018), respectively. Thus, we reasonably put the numerical wave generation and absorption at a location 15 and 6 m from the toe of fore-reef slope, which is also the location of left boundary. Behind the reef flat, transmitted waves are allowed to run up on the back-reef beach, but they cannot overtop out of the computational domain due to a solid wall condition at the right boundary. In addition, we set the top boundary free to the atmosphere. For the two side faces, we employed the “empty” boundary in OpenFOAM to simulate the 2-D reef configurations. When solitary waves interact with the investigated laboratory reefs, strong turbulence is expected to be generated inside the domain where the wave breaks near the reef crest and propagates on the reef flat as a moving bore; thus we do not set the inflow boundary condition with desired turbulence characteristics for the LES at the wave generation boundary. Meanwhile, since both the laboratory reef surfaces are very smooth, the flow structure near the bottom is not resolved in our simulations, and we only impose a no-slip boundary condition at the reef surfaces by adjusting the velocity near the bottom to satisfy the logarithmic law of the wall.

Figure 2 Numerical grids and boundary conditions of the numerical domain.

Structured mesh is used to discretize the computational domain. The discretization is kept constant in the spanwise (y) direction (one layer of 20 and 10 mm for the Roeber, 2010, and Yao et al., 2018, scenarios, respectively) while it varies in the streamwise (x) direction to reduce the number of the total cells. From the left boundary to the toe of the fore-reef slopes, Δx decreases gradually from 100 to 20 mm for the Roeber (2010) scenario and from 24 to 8 mm for the Yao et al. (2018) scenario (see, e.g., Fig. 2b and c). The core region (see, e.g., Fig. 2d), covering from the fore-reef slope to the back-reef wall or beach, maintains constant cell sizes of Δx=20 mm and 8 mm for the two experiments, respectively. Grid refinement near the free surface (e.g., Fig. 2b and c) is conducted at the core region in the x direction by reducing the grid sizes to one-quarter of their original values, e.g., Δx=5 and 2 mm. For the vertical (z) direction, the grid size is initially set to be Δz=20 and 8 mm across the domain for the Roeber (2010) and Yao et al. (2018) scenarios, respectively. Grid refinement near the free surface (e.g., Fig. 2b and c) is also conducted across the domain by reducing the grid sizes to Δz=5 and 2 mm. The total computational mesh consists of 4.87 and 1.18 million cells for the Roeber (2010) and Yao et al. (2018) scenarios, respectively. The simulation duration is appointed to be 80 and 30 s to guarantee the arrival of the reflected waves at the most offshore wave gauge in both experiments. The time step is automatically adjusted during computation for a constant Courant number of 0.25. Via parallel computing, it takes approximately 16 and 2 d for the Roeber (2010) and Yao et al. (2018) scenarios, respectively, on a cluster server with 44 CPUs (Intel Xeon, E5-2696, 2.2 GHz). No notable improvement of the results could be found with further refinement of the grid size.

For LES modeling solitary wave breaking over reefs, it is crucial to examine the Reynolds number (Re) at the incipient breaking point where strong turbulence is generated. It could be calculated by $Re=u_{\mathrm{b}}(H_{\mathrm{b}}+h_{\mathrm{b}})/\mathit{\nu }$ with $u_{\mathrm{b}}=c_{\mathrm{b}}{H}_{\mathrm{b}}/h_{\mathrm{b}}$ and $c_{\mathrm{b}}=\sqrt{g\left(H_{\mathrm{b}}+h_{\mathrm{b}}\right)}$, where H_b, h_b, u_b and c_b are wave height, water depth, streamwise velocity and wave celerity at the breaking point, respectively. Re is estimated for all tested scenarios by using H_b=H₀ and h_b=h_r (i.e., ignoring wave shoaling on the fore-reef slope and assuming wave breaking at the reef edge), and the values are also given in Table 1. Since the near-wall eddies are not resolved in this study, the total required grid number is independent of Re (Pope, 2000). The ideal grid size of the LES model should be down to the Kolmogorov scale, which is not feasible due to the limitation of computational resources. To test the convergence of grid size, we take the experiment with smaller wave flume (i.e., Scenario 3 in Table 1) which requires finer grid resolution as an example. We only examine the grid across the reef profile (the aforementioned core region) where the effect of grid size is supposed to be most influential. Both grid sizes (Δx and Δy) ranging from 8 mm down to 1 mm are tested. The results in terms of the dimensionless free-surface elevation and streamwise velocity associated with the leading solitary wave in the inner reef flat (G7) are compared in Fig. 3. Only differences of less than 2 % in terms of wave and flow could be observed with the grid size varying from 2 to 1 mm, indicating that our selection of grid size $\mathrm{\Delta }x=\mathrm{\Delta }z=\mathrm{2}$ mm is sufficient for the current simulations.

Figure 3 Variation of the maximum dimensionless free-surface elevation (${\mathit{\eta }}_{\max }/H_{\mathrm{0}}>\mathrm{1}$) and streamwise velocity ($u_{\max }/\sqrt{gh}>\mathrm{1}$) at G7 with the grid size (Δx and Δy) across the reef for case 3 in Table 1.

To evaluate the performance of the model, the model skill value is adopted and calculated by Wilmott (1981):

$$\begin{array}{}\text{(14)}& \mathrm{skill}=\mathrm{1}-{\displaystyle \frac{\sum {\left|Y_{\mathrm{model}}-Y_{\mathrm{obs}}\right|}^{\mathrm{2}}}{\sum {\left(\left|Y_{\mathrm{model}}-\stackrel{\mathrm{‾}}{Y_{\mathrm{obs}}}\right|+\left|Y_{\mathrm{obs}}-\stackrel{\mathrm{‾}}{Y_{\mathrm{obs}}}\right|\right)}^{\mathrm{2}}}},\end{array}$$

where Y_model is the predicted value and Y_obs is the measured value. The overline indicates that the average value is taken. The higher the skill number (close to 1), the better performance of the numerical model.

3.3 Comparison between numerical and experimental results

Figure 4 compares the computed and the measured cross-shore distribution of the free-surface elevations (η) at different stages (t) for Scenario 1, where η is normalized by the offshore still water depth (h₀) and t is normalized by $\sqrt{h_{\mathrm{0}}/g}$. The incident solitary wave gets steeper on the fore-reef slope at $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{62.3}$ due to the shoaling effect. Then its front becomes vertical prior to breaking at $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{64.3}$. At $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{65.8}$, a plunging breaker occurs with air entrainment and splash-up near the reef edge. After that, the breaking wave starts to travel on the reef flat in the form of a propagating turbulent bore at $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{67.1}$. The bore shows a gradual reduction in amplitude and continues to propagate downstream on the reef flat at $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{76.3}$. The numerical results generally agree well with the laboratory measurements at all stages with the skill values larger than 0.85, indicating the robustness of the adopted model to address the solitary wave processes across the laboratory reef profile in the large wave flume. When comparing the predictions between our Navier–Stokes-equation-based model and a Boussinesq model adopted by Roeber (2010), it seems that our model better captures the steep slope of the near-breaking wave ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{64.3}$) and breaking wave ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{65.8}$).

Figure 4 Dimensionless free-surface elevations (η∕h₀) across the reef at different stages ($t/\sqrt{h_{\mathrm{0}}/g})$ for Scenario 1. Red lines – present simulations; blue lines – simulations from Roeber (2010); open circles – measurements from Roeber (2010); skill values are for the present simulations.

Figure 5 Time series of dimensionless free-surface elevations (η∕h₀) at different cross-shore distances from the wavemaker (D) for Scenario 1. Red lines – present simulations; blue lines – simulations from Roeber (2010); open circles – measurements from Roeber (2010); skill values are for the present simulations.

Figure 5 illustrates the computed and measured time series of dimensionless free-surface elevations (η∕h₀) at different cross-shore locations (D) for Scenario 1. It appears that the model reasonably simulates the transformation processes of the solitary wave on the fore-reef slope (D=35.9 and 44.3 m) and near the reef edge (D=50.4 m) with the skill values larger than 0.9. The skill values relatively decrease right after the incipient wave breaking point (D=57.9 m) and at the central reef flat (D=65.2 m). Such discrepancies may be primarily due to the air entrainment in measuring both the breaking wave and the moving bore (Roeber, 2010) as well as the air bubble effect on free-surface tracking by the VOF method. In addition, the second peaks in the time series are due to wave reflection from the back-reef wall, which are well predicted by the present model. Meanwhile, no notable difference could be found in view of the time-series predictions between the present model and the model of Roeber (2010), except at D=65.2 m where the bore amplitude decays in our simulation compared to that at D=57.9 m.

Figure 6 Time series of dimensionless streamwise velocity ($u/\sqrt{gh}$) at different cross-shore distances from the wavemaker (D) for Scenario 1. Red lines – present simulations; blue lines – simulations from Roeber (2010); open circles – measurements from Roeber (2010); skill values are for the present simulations.

Figure 7 Dimensionless free-surface elevations (η∕h₀) across the reef at different stages ($t/\sqrt{h_{\mathrm{0}}/g}$) for Scenario 2. Red lines – present simulations; blue lines – simulations from Roeber (2010); open circles – measurements from Roeber (2010); skill values are for the present simulations.

Figure 6 depicts the time series of streamwise velocity (u) at five cross-shore locations (D) for Scenario 1, in which u is normalized by the local shallow water wave speed ($\sqrt{gh}$). The model satisfactorily captures the measured velocity offshore (D=17.8 m), on the fore-reef slope (D=47.4 m), on the central reef flat (D=72.6 m) and near the shoreline (D=80.2 m). A transition from the subcritical flow ($u/\sqrt{gh}<\mathrm{1}$) to supercritical flow ($u/\sqrt{gh}>\mathrm{1}$) could be observed right after wave breaking (D=61.6 m), and a less satisfactory prediction (skill values = 0.76) at this location is probably again due to both the effect of air bubbles on both flow measurements in the experiments and free-surface tracking in the simulations. Overall, the adopted model outperforms the Boussinesq model of Roeber (2010) in view of the velocity predictions, particularly both near the breaking point (D=61.6 m) and the shoreline on the reef flat (D=80.2 m).

Figure 8 Time series of dimensionless free-surface elevations (η∕h₀) at different cross-shore distances from the wavemaker (D) for Scenario 2. Red lines – present simulations; blue lines – simulations from Roeber (2010); open circles – measurements from Roeber (2010); skill values are for the present simulations.

As previously introduced, the reef profile of Scenario 2 is identical to that of Scenario 1 except for a reef crest located at the reef edge. The cross-shore distribution of dimensionless free-surface elevations (η∕h₀) at different stages ($t/\sqrt{h_{\mathrm{0}}/g}$) for Scenario 2 is demonstrated in Fig. 7. The steepened shoaling wave on the fore-reef slope appears at $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{65.0}$ and its front becomes almost vertical prior to breaking at $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{66.5}$. The breaking wave begins to overtop over the reef crest ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{69.1}$), and it then collapses on the lee side of reef crest, resulting in a moving turbulent bore ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{72.5}$). The bore travels shoreward on the reef flat with the continuous damping of its magnitude ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{80.5}$). The skill values for all sampling locations in this scenario are larger than 0.9, implying that the adopted model is able to well address the solitary wave processes over a more complicated reef geometry such as the presence of a reef crest at the reef edge. Again, the present model predicts the near-breaking wave ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{66.5}$) and breaking wave ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{69.1}$ and $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{72.5}$) slightly better than the model adopted by Roeber (2010).

Figure 8 compares the measured and simulated times-series of dimensionless free-surface elevations (η∕h₀) at various cross-shore locations (D) for Scenario 2. The skill values at all locations are larger than 0.85. It suggest again that the present model not only reasonably reproduces wave propagation offshore (D=17.6 m), shoaling on the fore-reef slope (D=35.9 and 44.3 m) and near breaking in front of the reef crest (D=50.4 m), breaking-wave transformation over the reef crest (D=57.9 m) and bore evolution on the reef flat (D=65.2 m), but also captures the tail waves caused by wave reflection from the back-reef wall (see, e.g., D=65.2 m). Overall, both our model and the model of Roeber (2010) reproduce the time series of free-surface elevations equally well for this scenario.

Figure 9 Time series of dimensionless streamwise velocity ($u/\sqrt{gh}$) at the cross-shore distance D=54.4 m from the wavemaker for Scenario 2. Red lines – present simulations; blue lines – simulations from Roeber (2010); open circles – measurements from Roeber (2010); skill values are for the present simulations.

As for Scenario 2, Roeber (2010) only reported one location of flow measurement on the seaside face of the reef crest. Figure 9 presents the time series of dimensionless streamwise velocity ($u/\sqrt{gh}$) at the point x=54.4 m, and a skewed and peaky velocity profile is observed to be associated with the leading solitary wave because the position is very close to the incipient wave breaking point. The two secondary peaks in the time series are generated by the reflected waves from the reef crest and from the back-reef wall, respectively. The model captures the temporal variation of the current fairly well with the skill value of 0.86, and its prediction is also better than that from the model of Roeber (2010), particularly for the reflected waves.

Figure 10 Time series of dimensionless free-surface elevations (η∕h₀) at different cross-shore sampling locations (G1–G8) for Scenario 3. Red lines – present simulations; blue lines – simulations from Yao et al. (2008); black lines – measurements from Yao et al. (2008); skill values are for the present simulations.

Figure 11 Time series of dimensionless free-surface elevations (η∕h₀) at different cross-shore sampling locations (G1–G8) for Scenario 4. Red lines – present simulations; blue lines – simulations from Yao et al. (2008); black lines – measurements from Yao et al. (2008); skill values are for the present simulations.

The experiments of Yao et al. (2018) only measured the time series of wave records at limited locations (G1–G8) across the reef as well as the maximum wave run-up on the final beach. Figure 10 compares the computed and measured time series of free-surface elevations for Scenario 3. The overall agreement between the simulations and experiments for G1–G8 is very good, with the skill values at all locations larger than 0.9. When the solitary wave travels from the toe (G2) to the middle of fore-reef slope (G3), it gets steeper due to the shoaling effect. Wave breaking starts at a location right before the reef edge (G4), and the surf zone processes extend over the reef flat in the form of a moving bore. Thus, from G5 to G8, the wave time series show saw-shaped profiles and there is a cross-shore decrease in the leading solitary wave height. Such features of the breaking waves are also well captured by the model. Note that the second peak in the time series of G7 is due to wave reflection from the back-reef beach, and the incident and reflected waves are not fully separated from each other at G8 because this location is too close to the beach. The predicted and measured wave run-ups are 0.122 and 0.109 m, respectively, for this scenario. Compared to the Boussinesq model employed by Yao et al. (2018), no significant difference in the predicted time series could be found for the present Navier–Stokes-equation-based model.

Figure 11 depicts the same comparison of wave time series but for the reef profile with a lagoon (Scenario 4). Again, the model performance for this scenario is fairly good (all skill values larger than 0.9). The predicted and measured wave run-ups are 0.123 and 0.116 m, respectively, for this scenario. Notable mismatch only appears for those small wave oscillations generated by the reflected wave propagating out of the lagoon to the reef flat (i.e., from G8 to G6). But our model seems to be superior to the model of Yao et al. (2018) in reproducing those oscillations at G7 and G8. We finally remark that the tail of the leading solitary wave, particularly from G1 to G4, is below the initial water level in the laboratory data, which is due to the water lost to form the generated wave crest around the paddle of the wavemaker. However, such a phenomenon is not observed in the numerical results because we generate a theoretical solitary wave in the numerical domain as indicated by Eq. (11).

4 Model applications

4.1 Effects of reef morphology variations on the solitary wave run-up

In this section, we apply the well-validated LES model to examine the variations of reef morphological parameters (fore-reef slope, back-reef slope, lagoon width, reef-crest width) that may affect the wave run-up (R) on the back-reef beach. Based on Scenario 3 (1:6 for both the slopes of fore reef and back reef, 9.6 m for the reef length, no reef crest and no lagoon) from Yao et al. (2018), we firstly test five slopes (1:2, 1:4, 1:6, 1:8 and 1:10, which all fall within the common range of 1:1 to 1:20 in the reported field observations; see, e.g., Table 1 of Quataert et al., 2015) for both the fore reef and the back reef. We then consider the existence of a lagoon at the rear of reef flat by testing four upper widths of the lagoon (1.6, 3.2, 4.8 and 6.4 m) and comparing to the case without a lagoon (lagoon width = 0 m). We finally investigate a trapezoidal reef crest located at the reef edge with its seaward slope matching the fore-reef slope and its shoreward slope of 1:1. We examine five reef-crest widths (0.1, 0.2, 0.3, 0.4 and 0.5 m) in view that the dimension of reef crest at the field scale is on the magnitude of meters (see, e.g., Hench et al., 2008). During simulations, each run is performed by changing one of above morphological parameters while keeping other parameters unaltered. All runs are conducted under a combination of one solitary wave height (H₀=0.08 m) and two reef-flat water depths (h_r=0.05 m and h_r=0.1 m).

Figure 12 The predicted wave run-up on the back-reef beach (R) under H₀=0.08 m with varying (a) fore-reef slopes, (b) back-reef slopes, (c) lagoon widths and (d) reef-crest widths.

Figure 13 Comparison of wave-driven current and vorticity on the x–z plane at different stages ($t/\sqrt{h_{\mathrm{0}}/g}$) between the reefs with and without the reef crest (H₀=0.08 m and h_r=0.05 m).

Generally, Fig. 12a shows that R is not very sensitive to the change in the fore-reef slope within the tested range, in that wave breaking for this scenario occurs close to the reef edge (G4); thus most of the surf zone processes and associated energy dispassion complete on the reef flat. Only when the fore-reef slope becomes steeper than 1:8, R decreases slightly under both water depths (h_r), which is attributed to the increased fore-reef reflection of the incident wave energy. Figure 12b reveals that R is more sensitive to the back-beach slope under both h_r. It decreases significantly as the back-beach slope becomes milder, which is consistent with that found for the plane slope (see, e.g., Synolakis, 1987). Figure 12c shows the variation of R with the lagoon width. Note that the zero width represents the reef without a lagoon. It appears that R increases notably with the increase in lagoon width because a wider lagoon dissipates less wave energy partly due to the stoppage of propagating bore and partly due to the reduction of bottom friction. As for the effect of reef-crest width (Fig. 12d), although the presence of a reef crest is reported to be an important factor affecting the wind wave transformation over fringing reefs (e.g., Yao et al., 2017), it seems to have little impact on the solitary wave run-up under both h_r; a slight decline of R could only be found under crest widths larger than 0.4. This is because the solitary wave is very long compared to the reef-crest width; thus most of its energy could transmit over the narrow reef crest. However, when the reef crest becomes sufficiently wide, its shallower crest tends to energize the wave breaking and thus the energy dissipation. To summarize all above analyses, it can be concluded that coastal regions protected by the fringing reefs with steeper back-reef slopes and wider lagoons are more valuable to coastal inundation during a tsunami event.

4.2 Wave-driven current and vortices around the reef crest and the lagoon

One advantage of the current Navier–Stokes-equation-based model over the depth-integrated models is its ability to resolve the vertical flow structure under breaking waves, particularly around the complex reef geometry. Based on the reef profile of Yao et al. (2018), Fig. 13 shows the simulated wave-driven current and vorticity on the x–z plane at different stages ($t/\sqrt{h_{\mathrm{0}}/g}$) for the reefs with and without a reef crest at the reef edge subjected to the same solitary wave condition (H₀=0.08 m and h_r=0.05 m). Without the reef crest, the shoaling wave propagates onto the horizontal reef flat with a uniform velocity distribution underneath ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{25.9}$ and 26.9), which is typical for the shallow water long waves. Until $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{27.9}$, wave breaking occurs in the form of a plunging breaker, and vortex generation gathers mainly around the wave crest. The vortices are transported further downstream at $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{28.9}$. When the wave crest exists, incipient wave breaking moves seaward and it takes place at the seaside edge of the reef crest ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{25.9}$). The breaker then overtops over the reef crest ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{26.9}$) and plunges onto the reef flat lee side of the reef crest, resulting in a hydraulic jump ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{27.9}$). Consequently, the wave-driven current at the rear part of the reef crest is dramatically increased compared to the same location without the crest. Both the intensity and the extent of vortex generation are also enlarged at the lee side of the reef crest ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{28.9}$), leading to increased wave energy dissipation compared to the case without the reef crest.

Figure 14 Comparison of wave-driven current and vorticity on the x–z plane at different stages ($t/\sqrt{h_{\mathrm{0}}/g}$) between the reefs with and without the lagoon (H₀=0.08 m and h_r=0.05 m).

Figure 14 compares the computed wave-driven current and vorticity on the x–z plane at different stages ($t/\sqrt{h_{\mathrm{0}}/g}$) between the reefs in the presence and absence of the lagoon. Without the lagoon, the propagating bore arrives with strong vortex motions ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{49.4}$). The vortices are eventfully transported downstream from $t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{54.4}$ to 64.4. However, when the lagoon is present, the current speed over the depth slows down and additional vortices generate at the seaside edge of the lagoon as the bore propagates into the lagoon ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{49.4}$). The peak value of the vorticity appears at a later time ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{54.4}$). After that, the vortices in the lagoon are primarily diffused by the vortex shedding ($t/\sqrt{h_{\mathrm{0}}/g}=\mathrm{59.4}$ and 64.4). Compared to the case without the lagoon, although the existence of a lagoon dissipates less wave energy by terminating the propagating bore and reducing the reef-flat friction as previously stated, the vortex generation and diffusion in the lagoon as demonstrated here are believed to cause local energy loss. We finally remark that the wave-driven current and vortices examined in this section could provide a first step to analyze more sophisticated problems, such as the tsunami-induced sediment/debris transport over the fringing reefs.

5 Conclusions

To remedy the inadequacies of using the depth-integrated models to simulate the interaction between tsunami-like solitary waves and fringing reefs, a 3-D numerical wave tank, solving the Navier–Stokes equations with the LES for turbulence closure, has been developed based on the open-source CFD tool OpenFOAM^®. The free surface is tracked by the VOF method. Two existing laboratory experiments with the wave, flow and wave run-up measurements based on different fringing reef profiles are employed to validate the numerical model. Simulations show that the current Navier–Stokes-equation-based model outperforms the commonly used Boussinesq-type models in view of its capability to better reproduce the breaking waves and wave-driven current on the reef flat. The model is then applied to investigate the impacts of varying morphologic features on the back-reef wave run-up. The flow and vorticity fields associated with the breaking solitary wave around the reef crest and the lagoon are also analyzed via the numerical simulations.

Model results shows that wave run-up on the back-reef slope is most sensitive to the variation of the back-reef slope, less sensitive to the lagoon width, and almost insensitive to the variations of both the fore-reef slope and the reef-crest width within our tested ranges. The existence of a reef crest or a lagoon can notably alter the wave-driven current and vortex evolutions on the reef flat. These findings demonstrate that low-lying coastal areas fringed by coral reefs with steep back-reef slopes and larger lagoons are expected to experience larger wave run-up near the shoreline; thus they are more susceptible to coastal inundation during a tsunami event.