Numerical Investigations of Tsunami Run-Up and Flow Structure on Coastal Vegetated Beaches

Abstract

Tsunami waves become hazardous when they reach the coast. In South and Southeast Asian countries, coastal forest is widely utilized as a natural approach to mitigate tsunami damage. In this study, a depth-integrated numerical model was established to simulate wave propagation in a coastal region with and without forest cover. This numerical model was based on a finite volume Roe-type scheme, and was developed to solve the governing equations with the option of treating either a wet or dry wave front boundary. The governing equations were modified by adding a drag force term caused by vegetation. First, the model was validated for the case of solitary wave (breaking and non-breaking) run-up and run-down on a sloping beach, and long periodic wave propagation was investigated on a partially vegetated beach. The simulated results agree well with the measured data. Further, tsunami wave propagation on an actual-scale slope covered by coastal forest Pandanus odoratissimus (P. odoratissimus) and Casuarina equisetifolia (C. equisetifolia) was simulated to elucidate the influence of vegetation on tsunami mitigation with a different forest open gap. The numerical results revealed that coastal vegetation on sloping beach has significant potential to mitigate the impacts from tsunami waves by acting as a buffer zone. Coastal vegetation with open gaps causes the peak flow velocity at the exit of the gap to increase, and reduces the peak flow velocity behind the forest. Compared to a forest with open gaps in a linear arrangement, specific arrangements of gaps in the forest can increase the energy attenuation from tsunami wave. The results also showed that different cost-effective natural strategies in varying forest parameters including vegetation collocations, densities, and growth stages had significant impacts in reducing the severity of tsunami damage.

1. Introduction

Tsunamis are generated by marine earthquakes, underwater volcanic eruptions, or submerged landsides. Recent tsunami disasters (e.g., in the Indian Ocean in 2004, Samoa in 2009, Chile in 2010, Tohoku in 2011, and Indonesia in 2018) have caused vast losses of life and property, destruction of critical infrastructure in low-lying coastal areas, and massive damage to the coastal ecosystem [1,2,3,4]. These great threats demonstrate the need to mitigate tsunamis by using artificial obstacles including wave-dissipating blocks, rock breakwaters, and large embankments. Recent studies investigated tsunami waves in both laboratory experiments [5,6,7,8] and field investigations [9,10] with the ultimate goal of mitigating their impacts and reducing natural disasters in coastal regions. Considering the advantages of flexibility and repeatability in numerical models, many numerical simulations have been carried out using different models to investigate the run-up heights and velocity processes of tsunamis [11,12,13,14,15,16]. Typical models that were considered included depth-integrated Boussinesq-type models and depth-integrated non-linear shallow water models (NSWMs), with results showing that NSWMs are more robust and efficient for tsunami predictions in real-scale basins.

Generally, coastal forests which act as natural obstacles are widely distributed in coastal areas. They also effectively mitigate tsunami damage from economic, environmental, and aesthetically pleasing points of view in concert with hard structures [17,18]. Currently, there is some literature on the capability of coastal vegetation to attenuate wave energy in short-period waves [19,20]. While tsunamis are categorically different from short-period waves, there have been an increasing number of studies that address the role of vegetation in mitigating coastal natural disasters, which can include the strategy of planting coastal vegetation as a bio-shield [21,22,23]. Recent research has calculated the friction and drag force of vegetation in a tsunami, using mathematical equations. Real-scale simulations of tsunami were performed to investigate the effectiveness of forests as a bio-shield for tsunami protection, with results indicating that a forest with two layers in the vertical direction including P. odoratissimus and C. equisetifolia was effective for attenuating the wave energy [24]. Based on two dimensional nonlinear long-wave equations, Thuy et al. developed a numerical model to estimate tree-breaking, the drag forces, and turbulence-induced shear force due to the presence of vegetation (P. odoratissimus and C. equisetifolia) [25]. Considering the effects of the porosity of vegetation, a one-dimensional numerical model using Boussinesq-type equations was developed to evaluate the effects of the forest density distribution on tsunami-force reduction, including the drag and inertia forces caused by vegetation [26]. Numerical investigations were conducted to study the impact of patchy vegetation on tsunami dynamics based on Boussinesq model, and the results demonstrated that patchy vegetation, with appropriate configuration, can be effective in mitigating tsunami hazards [27]. Based on a nonlinear long wave equation model that included the breaking or washout of trees, numerical simulations were carried out to estimate the effects of coastal forest and sea embankments on reducing the washout area and the tsunami mitigation function of the coastal forest [28].

Roads are perpendicular to the coast function as pathways that form gaps in the vegetation of the coastal forest. Gaps in the coastal forest can increase risks and potential damage, as the water flow from the tsunami accelerates as it moves through the gap into the densely populated block [29]. Tanimoto et al. studied it by changing the width of open gap and found that a specific gap with 15 m width causes the highest flow velocity in their simulated conditions [30]. Thuy et al. numerically simulated the effect of an open gap in coastal forests on tsunami run-up, and found that maximum flow velocity greatly increased at the open gap exit, meaning that an open gap (like a road) in a coastal forest had a negative effect on tsunami run-up behind the forest [31]. Nandasena et al. investigated hydrodynamic parameters (flow velocity and depth) of tsunami waves and bending moment of vegetation on Misawa, a site covered by pine forest with two gaps [32]. However, no studies have discussed the effects of different open gaps in existing forests via the analysis of maximum tsunami run-up and variations of velocity through the gaps using a model with shock capture capability.

In this paper, we numerically investigate the hydrodynamic processes on a vegetated sloping beach and quantify the effects of a gap in forest vegetation on the run-up by solving the depth-averaged 2D model. The proposed model was tested for breaking and non-breaking solitary waves propagating on a bare sloping beach, and long periodic wave propagation on a partially-vegetated sloping beach to examine the accuracy of the numerical model. Then, the model was used to simulate tsunami waves propagating on actual-scale forest-covered beach, and to study the mitigation effects of gap arrangements and different vegetation parameters on tsunami waves.

2. Numerical Method

2.1. Governing Equations

The depth-averaged 2D shallow water equations are formed by integrating the Navier–Stokes equations which include continuity and momentum equations for depth-averaged free surface flows. These are expressed as Equations (1)–(3):

$$\frac{\partial h}{\partial t}+\frac{\partial uh}{\partial x}+\frac{\partial vh}{\partial y}=0$$

(1)

$$\frac{\partial hu}{\partial t}+\frac{\partial (huu)}{\partial x}+\frac{\partial (hu\mathrm{v})}{\partial y}-\frac{\partial }{\partial x}(v_{t}h\frac{\partial u}{\partial x})-\frac{\partial }{\partial y}(v_{t}h\frac{\partial u}{\partial y})=-gh\frac{\partial \eta }{\partial x}-{\tau }_{bx}+f_{c}hv-f_x$$

(2)

$$\frac{\partial hv}{\partial t}+\frac{\partial (huv)}{\partial x}+\frac{\partial (hvv)}{\partial y}-\frac{\partial }{\partial x}(v_{t}h\frac{\partial v}{\partial x})-\frac{\partial }{\partial y}(v_{t}h\frac{\partial v}{\partial y})=-gh\frac{\partial \eta }{\partial y}-{\tau }_{by}-f_{c}hu-f_y$$

(3)

In the above equations, t is time, h indicates the local water depth, u and v stand for the depth-averaged flow velocities in the x and y directions, respectively, η stands for the water surface elevation from a reference datum, v_t is defined as the eddy viscosity coefficient calculated by ${\nu }_t=\alpha u_{\ast }h$ where α is empirical constant and ranges from 0.3 to 1.0, $u_{\ast }$ means the bed shear velocity, ${\tau }_{bx}$ and ${\tau }_{by}$ denote friction in the x and y directions, respectively, ${\tau }_{bx}=\mathrm{g}\frac{n^{2}u\sqrt{u^2+v^2}}{h^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$, ${\tau }_{by}=\mathrm{g}\frac{n^{2}v\sqrt{u^2+v^2}}{h^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$, where n stands for Manning’s roughness coefficient, f_c indicates the Coriolis parameter, and f_x and f_y stand for the drag force caused by vegetation.

The conservative and vector form of 2D shallow water equations with vegetation effects are expressed as follows:

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x}+\frac{\partial \mathbf{G}}{\partial y}=\frac{\partial {\mathbf{F}}_d}{\partial x}+\frac{\partial {\mathbf{G}}_d}{\partial y}+\mathbf{S}$$

(4)

The definitions of U, F, G, F_d and G_d are shown as follows:

$$\mathbf{U}=\left[\begin{array}{c}h\\ hu\\ hv\end{array}\right],\mathbf{F}=\left[\begin{array}{c}hu\\ huu\\ huv\end{array}\right],\mathbf{G}=\left[\begin{array}{c}hv\\ huv\\ hvv\end{array}\right],{\mathbf{F}}_d=\left[\begin{array}{c}0\\ {\nu }_t\frac{\partial uh}{\partial x}\\ {\nu }_t\frac{\partial vh}{\partial x}\end{array}\right],{\mathbf{G}}_d=\left[\begin{array}{c}0\\ {\nu }_t\frac{\partial uh}{\partial y}\\ {\nu }_t\frac{\partial vh}{\partial y}\end{array}\right],\mathbf{S}=\left[\begin{array}{c}0\\ -gh\frac{\partial \eta }{\partial x}-{\tau }_{bx}+f_{c}hv-f_x\\ -gh{\frac{\partial \eta }{\partial y}}_{}-{\tau }_{by}-f_{c}hu-f_y\end{array}\right]$$

(5)

For notational convenience, Equation (4) is often rewritten as:

$$\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot {\mathbf{E}}_w^{adv}=\nabla \cdot {\mathbf{E}}_w^{dif}+\mathbf{S}$$

(6)

where ${\mathbf{E}}_w^{adv}=\mathbf{F}i+\mathbf{G}j$ and ${\mathbf{E}}_w^{dif}={\mathbf{F}}_{d}i+{\mathbf{G}}_{d}j$.

2.2. Vegetation Drag Force

The resistance effects of vegetation on flow are added into the momentum equations as an internal source of resistant force per unit fluid mass. The drag force exerted by vegetation per unit volume is expressed as [33]:

$$f_x=\frac{1}{2}N{C}_D\left(h\right)bv\min (h_v,h)u\sqrt{u^2+v^2},\hspace{1em}{f}_y=\frac{1}{2}N{C}_D\left(\mathrm{h}\right)b_v\min (h_v,h)v\sqrt{u^2+v^2}$$

(7)

Here, N stands for the vegetation density defined as the number of plants per square meter, b_v indicates the diameter, and h_v indicates the height of the vegetation. C_D (h) indicates the depth-averaged equivalent drag coefficient which considers the vertical stand structures of a tree, defined by Tanaka et al. [1] as follows:

$$C_D\left(h\right)=C_{D-ref}\frac{1}{h}{\displaystyle \underset{0}{\overset{\mathrm{h}}{\int }}\alpha \left(z_g\right)}\beta \left(z_g\right)d{z}_g$$

(8)

$$\alpha \left(z_g\right)=\frac{b\left(z_g\right)}{b_{ref}}$$

(9)

$$\beta \left(z_g\right)=\frac{C_D\left(z_g\right)}{C_{D-ref}}$$

(10)

where C_D-ref is the reference drag coefficient of the trunk at z_g (equal to 1.2 m in principle), α(z_g) is considered as an additional coefficient to express the effects of cumulative width on drag force at each height z_g, b(z_g) means the projected width, and b_ref indicates the reference projected width. β(z_g) is expressed as an additional coefficient representing the effect of leaves or aerial roots on drag force, and C_D(z_g) is drag coefficient of a tree at the height z_g above the ground surface.

2.3. Finite Volume Method

The discretization of the governing equations is based on the finite volume method using an unstructured triangular mesh. The conserved variables are defined at the cell centers and represent the average value of each cell. The integral form of Equation (10) over the ith control volume can be expressed as:

$${\displaystyle \underset{V_i}{\int }\frac{\partial \mathbf{U}}{\partial t}}{d}_V+{\displaystyle \underset{V_i}{\int }\nabla \cdot {\mathbf{E}}_w^{adv}}{d}_V={\displaystyle \underset{V_i}{\int }\nabla \cdot {\mathbf{E}}_w^{dif}}{d}_V+{\displaystyle \underset{V_i}{\int }\mathbf{S}}{d}_V$$

(11)

In this expression, subscript i and j denote the element and the element side, respectively, while V_i stands for the domain of the ith. Basing Green’s theorem. Equation (11) can be rewritten as:

$$\frac{\Delta {\mathbf{U}}_i}{\Delta t}{A}_i=-{\displaystyle \underset{L_i}{\oint }{\mathbf{E}}_w^{adv}\cdot \mathbf{n}}dl+{\displaystyle \underset{L_i}{\oint }{\mathbf{E}}_w^{dif}\cdot \mathbf{n}}dl+{\displaystyle \underset{V_i}{\int }\mathbf{S}}{d}_V$$

(12)

where ${\mathrm{U}}_i$ stands for the average value of the conserved variables over the ith cell, and ${\mathrm{U}}_i$ is stored at the center of the ith cell, with ${\mathbf{U}}_i=\frac{1}{A_i}{\displaystyle \underset{V_i}{\int }\mathbf{U}}{d}_V$. L_i means the boundary of the V_i. $A_i$ denotes the area of the ith cell, n is the outward surface normal vector of L_i, with $\mathbf{n}=(n_x,n_y)=(\cos \varphi ,\sin \varphi )$, and $\varphi$ is the angle included between the x direction and the outward normal vector.

In the 2D triangular grid system, the line integral term in the Equation (12) can be further approximated and assessed as follows:

$$\Delta {\mathbf{U}}_i=-\frac{\Delta t}{A_i}{\displaystyle \sum _{j=1}^m\left({\mathbf{E}}_w^{adv}\cdot {\mathbf{n}}_{ij}\right)l_{ij}+\frac{\Delta t}{A_i}{\displaystyle \sum _{j=1}^m\left({\mathbf{E}}_w^{dif}\cdot {\mathbf{n}}_{ij}\right)}}{l}_{ij}+\frac{\Delta t}{A_i}{\displaystyle \underset{V_i}{\int }\mathbf{S}}{d}_V$$

(13)

where m denotes the number of total edges of the triangular cell (three in this model), the subscript j means the index of the edge of a triangular mesh, ${\mathbf{n}}_{ij}$ stands for the outward normal flux vector, and $l_{ij}$ means the length of the arc.

2.4. Evaluation of Numerical Fluxes

Variables are usually approximated as constant states within each control volume. Riemann problems at the interface of the cell can be solved using various Riemann approximations for evaluating the interface fluxes. The interface fluxes of Roe’s solver are expressed as follows:

$${\mathbf{E}}_w^{adv}\cdot \mathbf{n}=\frac{1}{2}\left[{(\mathbf{F},\mathbf{G})}_R\cdot \mathbf{n}+{(\mathbf{F},\mathbf{G})}_L\cdot \mathbf{n}-\left|J\right|({\mathbf{U}}_R-{\mathbf{U}}_L)\right]$$

(14)

where ${\mathbf{U}}_R$ and ${\mathbf{U}}_L$ are reconstructed Riemann state variables on the right and left sides of the cell interface, respectively. The flux Jacobian matrix A can be assessed as:

$$\left|J\right|=\frac{\partial \left({\mathbf{E}}_w^{adv}\cdot \mathbf{n}\right)}{\partial \mathbf{U}}=\frac{\partial \mathbf{F}}{\partial \mathbf{U}}{n}_x+\frac{\partial \mathbf{G}}{\partial \mathbf{U}}{n}_y=\left[\begin{array}{ccc}0& n_x& n_y\\ \left(c^2-u^2\right)n_x-uv{n}_y& 2u{n}_x+v{n}_y& u{n}_y\\ -uv{n}_x+\left(c^2-v^2\right)n_y& v{n}_x& u{n}_x+2v{n}_y\end{array}\right]$$

(15)

where $n_x$ and $n_y$ denote the components of the outward surface normal vector in the x- and y- directions, respectively, and c indicates the wave velocity, with

$$c=\sqrt{gh}$$

where $\left|J\right|=R\left|\Lambda \right|L$, R and L stand for the right and left eigenvector matrices, and $\left|\mathsf{\Lambda }\right|$ denotes the diagonal matrix of the absolute values of the eigenvalues of A.

Where $\left|\mathsf{\Lambda }\right|$ is defined as:

$$\left|\mathsf{\Lambda }\right|=\left(\begin{array}{ccc}{\mathsf{\Lambda }}^1& 0& 0\\ 0& {\mathsf{\Lambda }}^2& 0\\ 0& 0& {\mathsf{\Lambda }}^3\end{array}\right)$$

where

$${\lambda }^1=u{n}_x+v{n}_y,{\lambda }^2=u{n}_x+v{n}_y-c,{\lambda }^3=u{n}_x+v{n}_y+c$$

(16)

The right and left eigenvector matrices are expressed as follows:

$$R=\left(\begin{array}{ccc}0& 1& 1\\ n_y& u-c{n}_x& u+c{n}_x\\ -n_x& v-c{n}_y& v-c{n}_y\end{array}\right),L=\left(\begin{array}{ccc}-(u{n}_y-v{n}_x)& n_y& -n_x\\ \frac{u{n}_x+v{n}_y}{2c}+\frac{1}{2}& \frac{-n_x}{2c}& \frac{-n_y}{2c}\\ -\frac{u{n}_x+v{n}_y}{2c}+\frac{1}{2}& \frac{n_x}{2c}& \frac{n_y}{2c}\end{array}\right)$$

(17)

where the Riemann state variables u, v, and $c$ on the cell interface are necessary to deal with the fluxes as calculated by Roe’s average:

$$u=\frac{\sqrt{h_{+}}{u}_{+}+\sqrt{h_{-}}{u}_{-}}{\sqrt{h_{+}}+\sqrt{h_{-}}},v=\frac{\sqrt{h_{+}}{v}_{+}+\sqrt{h_{-}}{v}_{-}}{\sqrt{h_{+}}+\sqrt{h_{-}}},c=\sqrt{\frac{g(h_{+}+h_{-})}{2}}$$

(18)

The subscripts $+$ and $-$ indicate the right and left sides of the cell edge, respectively. When the drying-wetting interfaces exist in the computational domain, Roe’s average can be calculated by: $u=\frac{u_{+}+u_{-}}{2}$, $v=\frac{v_{+}+v_{-}}{2}$.

2.5. Treatment of Wetting and Drying Fronts

A technique to treat both wet and dry boundaries was introduced to achieve zero mass error [34,35]. A criterion, ε, was adopted to define and classify the following four types of edges:

• Wet edge (see Figure 1a): two adjacent cells are wet, in which water depth of left cell h_L > ε and water depth of right cell h_R > ε.
• Partially wet edge (with flux), as presented in Figure 1b: a wet cell (left) links to a dry cell on the right, and the water level of the wet cell is higher than that of the dry cell, where h_L > ε, h_R ≤ ε and water level of left cell η_L > water level of left cell η_R.
• Partially wet edge (no flux), as shown in Figure 1c: a wet cell (left) links to a dry cell on the right, and the water level of the wet cell is lower than that of the dry cell, where h_L > ε, h_R ≤ ε, and η_L < η_R. To eliminate the non-physical flux problem produced in the interface, the water level η_R and bed level Z_bR for the dry cell were temporarily replaced by a value which equaled to the water level η_L in the wet cell.
• Dry edge, as presented in Figure 1d: two adjacent cells are dry, where h_L ≤ ε and h_R ≤ ε.

According to these four types of edges, cells were correspondingly divided into three types:

• Wet cell: all the edges of this cell consisted of a wet or partially wet edges (with flux) and all the nodes of the cell are flooded.
• Dry cell: all the edges of this cell consist of dry or partially wet edges (no flux).
• Partially wet cell: all other cells do not satisfy the criteria of either a wet or dry cell, as defined above.

3. Numerical Simulation and Experimental Validation

3.1. Solitary Wave Run-up on a Bare Sloping Beach

The depth-integrated shallow water model was run for experiment cases of breaking and non-breaking solitary waves on a bare sloping beach conducted by Synolakis to validate the accuracy of the numerical scheme for modeling wave run-up and run-down [36]. The topography consists of a 1:19.85 sloping beach adjacent to a constant depth region, as illustrated in Figure 2. A solitary wave propagates from left to right in a wave flume; and according to first-order solitary wave theory, the water surface elevation and velocity in initial time are defined as follows:

$$\eta \left(x\right)=H_W{\sec h}^2\left(\sqrt{\frac{3HW}{4{h_0}^3}}\left(x-X_0\right)\right)$$

(19)

$$u\left(x\right)=\eta \left(\mathrm{x}\right)\sqrt{\frac{g}{h_0}}$$

(20)

where H_w means wave height, h₀ denotes initial water depth with a value of 1 m, X₀ stands for the position of initial wave crest and is located a half wave length from the toe of the sloping beach in the computing domain, and u is wave velocity.

In this case study, the numerical model used uniform triangular cells with a grid spacing of 0.02 m and a time step of 0.001 s. The minimum water depth to define a “dry bed” was set as 0.0001 m, and the Manning’s bed roughness coefficient was calibrated as 0.01 for the glassed flume in the present model. The solitary wave was defined as non-breaking in this case, as H_w/h₀ = 0.0185. For convenience, the results are presented in non-dimensional forms: ${\mathrm{x}}^{\ast }=\frac{x}{h_0}$, ${\eta }^{\ast }=\frac{\eta }{h_0}$ and ${\mathrm{t}}^{\ast }=t\sqrt{\frac{\mathrm{g}}{h_0}}$. Comparisons of the simulated and experimental free-surface evolutions are presented in Figure 3. As illustrated, the incident wave propagates on the sloping beach at the early stage (t^∗ = 25, 30, 35, 40, 45, and 50), and reaches maximum run-up height at about t^∗ = 55, at which point backwash occurs. The maximum run-down happens at around t^∗ = 70. Simulated free surface profiles show good agreement with experimental data. Figure 4 shows the comparison of simulated and experimental water surface processes from a breaking solitary wave with H_w/h₀ = 0.3. In Figure 4, wave breaking is not well reproduced by the depth-averaged shallow model and the computed wave fronts are steeper and slightly earlier than the experimental results at t^∗ = 15 and t^∗ = 20, as the model does not consider a wave dispersion term [37,38]. However, the wave breaking is only in a small portion of the domain. At the next step (t^∗ = 25), wave breaking is reasonably simulated by the model as a collapse of the wave near the shoreline approximately which can be explained that the wave breaking may be weakened and become smaller portion of the domain. The breaking solitary wave reaches a maximum height around t^∗ = 45, and approaches the lowest position where a hydraulic jump is formed near the shoreline around t^∗ = 55. The simulated results agree with experimental data very well, which indicate that the proposed model can accurately predict the propagation of breaking and non-breaking solitary waves on a bare sloping beach.

3.2. Propagation of Long Periodic Waves on a Partially Vegetated Sloping Beach

A laboratory experiment exploring the propagation of long periodic waves in a wave channel with a partially-vegetated sloping beach was conducted in Saitama University [31,32]. In this experiment, the wave tank was 15 m long and 0.4 m wide for propagation, and the bed profile of the beach consisted of several sections with various slopes (Figure 5a). The vegetation was modelled using wooden cylinders with a diameter of 0.005 m at a still water depth of 0.44 m. The vegetation domain was from x = 10.36 to 11.36 m, with different vegetation gap widths (B_g = 0, 0.07, and 0.4 m), while the vegetation density was set as 2200 stem/m⁻², and the drag coefficient C_D of wooden cylinders was calibrated as 2.5 in the present study. An incident sinusoidal wave with a period of T = 20 s and wave height of 0.16 cm was propagated from a 0.52 m long constant bottom segment to a sloping beach. In the current case, water surface elevations were measured using capacitance wave gauges at six positions (G1–G6) along the center of the flume in cases of B_g = 0 and 0.4 m, and the flow velocities were measured using electromagnetic current meters at locations in the cross section passing behind vegetation domain G6 (see Figure 5b) in case of B_g = 0.07 m. The run-up height above still water surface was measured by tracing the moving of water front by eyes with a scale on the slope. The computational domain was represented by uniform triangular cells with a grid spacing of 0.01 m and a time step of 0.002 s for the numerical model. The minimum depth criterion for a dry and wet bed was considered as 0.001 m, and the Manning coefficient was set as 0.012.

A comparison of wave crests, heights, and troughs in Figure 6a,b demonstrates that for the case of non-vegetated areas, wave heights increased due to shoaling as the long waves approached shallow water, however, the increased amplitudes of wave heights were attenuated by the resistance effects of vegetation at the sloping zone, with the attenuated rate reaching 38.96% at G6 (at the center behind vegetation). The wave height at the front of the vegetated zone was also affected by the reflection of waves caused by the vegetation. Figure 7 shows the comparison of time series of velocities at the center of the gap exit and the center of the end of the vegetation when B_g = 0.07 m, indicating that the peak flow velocity at the center of the gap exit can reach 0.42 m/s and is 3.07 times more than the peak flow velocity at behind the center of the vegetation zone. Figure 8 shows the distribution of the peak velocities averaged from five wave periods at steady state in a cross-section of Gage 6, indicating that vegetation plays an important role in the distribution of flow velocities. The results demonstrate that the present model is an effective tool to predict long periodic wave propagation on a partially vegetated sloping beach. Vegetation can effectively attenuate the wave propagation; however, an open gap in vegetation zone generates large flow and adverse effect at gap exit.

3.3. Effects of Forest on Tsunami Run-up at Actual Scale

3.3.1. Coastal Topography and Forest Conditions

Coastal forests are widely considered to mitigate tsunami damage on coastal beaches in some countries and regions. In this case study, a numerical simulation was performed to investigate the run-up of tsunami and the effects of forest vegetation on tsunami mitigation. A coastal topography adopted for tsunami simulation is presented in Figure 9 [25,31]. The bed profile of the computed domain consisted of four slopes: S = 1/10, 1/100, 1/50, and 1/500, while the offshore water depth with a horizontal bed was 102 m below the still water level, all parameters for a tsunami wave simulation can be seen in references [25,31].

P. odoratissimus is a representative tree of South and Southeast Asia with a complex aerial root structure growing in coastal beach. In addition, C. equisetifolia is another representative species that grows densely on sandy beaches. These two species are discussed in evaluating the resistance of vegetation on tsunami propagation. The tree heights h_v of P. odoratissimus and C. equisetifolia are considered as 8 m and 11 m, respectively; the reference diameters b_ref of P. odoratissimus and C. equisetifolia are set as 0.2 m and 0.15 m, respectively; and the density values N of P. odoratissimus and C. equisetifolia are 0.22 trees/m² and 0.4 trees/m², respectively. The depth-averaged equivalent drag coefficient C_D(h) are variable with vegetation parameters and inundation depth h, and they are modified by Tanaka et al. and Thuy et al. [22,26]. The coastal forest starts at a slope of 1/500 (see Figure 9).

In this numerical simulation, the computational domain is discretized into 38,010 triangular meshes, in which the coarsest grid spacing at the sea-side is 20 m, and the finest grid spacing is 2.5 m at the nearshore area. The Manning’s roughness coefficient is set as 0.02 for the computational domain. By comparing the maximum water surface elevations in Figure 10, there is good agreement between the simulated and Thuy’s results [31].

In non-vegetated areas, as the wave crest approaches a 1/50 slope, its energy compresses and the maximum water surface elevations increase slightly because of shoaling. When the wave crest approaches a slope of 1/500, the maximum water surface elevations decrease due to the accelerated flow velocity on a mild slope. For fully vegetated areas, the maximum water surface elevations increase significantly at the front of the forest due to the reflection of the waves caused by the vegetation. In the forest zone, the maximum water surface elevations decrease drastically. The run-up height above the still water level is 5.10 m with full vegetation, which is lower than the 6.86 m for non-vegetation areas. The results indicate that P. odoratissimus effectively protects coastal communities from the destruction by a tsunami wave by reducing wave energy.

3.3.2. Effects of Forest with a Straight Open Gap on Tsunami Run-up

Figure 11 shows four types of gap arrangements (Cases 1, 2, 3, and 4) in P. odoratissimus forest with the simulation conditions of L_F = 200 m, B_F = 200 m, and B_G = 15 m. In this section, tsunami run-up on the sloping beach is considered using a straight gap (Case 1) to investigate the effect on tsunami run-up.

Figure 12a–c presents the predicted flow structure in Case 1 at different times (600, 660 and 700 s). The flow velocity of the first wave is fast in an open gap zone, while the wave moves slowly in the forest domain due to the resistance from vegetation. The first wave front quickly passes the open gap and spreads out from the exit, and the flow behind the forest gradually becomes uniform, with the exception of the areas close to the gap exit. Figure 13 shows the temporal variation of water surface elevations and flow velocities at four predicted stations (G1, G2, G3, and G4) for Case 1. In this case, the maximum water surface elevations in the presence of vegetation are higher at G1 (the middle of the open gap) and G2 (the middle before the vegetation) than those in the absence of vegetation due the reflection of tsunami waves caused by P. odoratissimus forest.

As the flow passes the forest, the maximum water surface elevations drastically decrease. The maximum water surface elevation values at G3 (the middle behind the forest) and G4 (the middle at the end of the open gap) in the presence of vegetation with an open gap decrease to 80.72% and 82.52% of that for the case of non-vegetation, respectively. The time delays in a tsunami arriving at G3 and G4 are about 45 and 14 s, respectively, when the wave passes through forest with an open gap compared to that of non-vegetated case. The P. odoratissimus forests slightly impact the peak velocities at G1 and G2, which reach 4.36 and 4.35 m/s, respectively. Due the resistance effects of P. odoratissimus, the flow velocities at G1 and G2 decrease dramatically after reaching their peak. The peak velocity at G3 is decreased by 35.13% and the peak velocity at G4 is increased by 44.04%. The open gap in forest vegetation amplifies negative effects on tsunami propagation even more than in the non-vegetated case.

3.3.3. Effects of Arrangements of Vegetation Patches and Open Gaps

Four types of gap structures are discussed to investigate their potential to increase the energy attenuation of tsunami waves on a vegetated beach. Figure 14 illustrates the structures of a flow field at 660 s for three cases with different gap arrangements, and shows that the flow velocities are larger in open gaps than those in vegetated zones, and the front of the first wave quickly passes the open gap, and the maximum peak flow velocities appear around the gap exit. The furthest distance that tsunami waves reach is variable, but there are differences in flow where structures occur. As presented in Figure 15, the effects of the gap arrangements on the water surface elevations at G3 and G4 are insignificant for four cases, while the different gap arrangements have important effects on flow velocities of tsunami waves.

Compared to Case 1, the peak flow velocities at the G3 station are increased by about 0.15 m/s (4.92%) and 0.21 m/s (6.71%) for Case 2 and Case 4, respectively, and the peak flow velocity is attenuated slightly by about 0.05 m/s (−1.46%) for Case 3. Compared to Case 1, the peak flow velocities at G4 are decreased by about 0.47 m/s (−6.77%), 0.34 m/s (−4.87%) and 0.89 m/s (−12.88%) for Case 2, 3, and 4, respectively. Figure 16 shows the comparisons of the run-up heights of the first tsunami wave for different cases including fully vegetated, non-vegetated, and four types of vegetation arrangements with different gaps. The results indicate that the slight differences of run-up heights occur for different gap arrangements, however, the presence of vegetation significantly attenuates the energy of tsunami waves compared to non-vegetated case. Therefore, it is reasonable to reduce the disadvantages from open gaps by changing gap arrangements.

3.3.4. Effects of Forest Parameters on Tsunami Run-up

To investigate the effects of coastal forest parameters on tsunami run-up, a sensitivity analysis of different tree configurations, densities, and growth stages was examined in the gap arrangement of Case 4. In this section, Collocation 1 (pure P. odoratissimus), Collocation 2 (pure C. equisetifolia), Collocation 3 (C. equisetifolia as the back vegetation layer of P. odoratissimus), and Collocation 4 (C. equisetifolia as the front vegetation layer of P. odoratissimus) are considered to analyze effects of tree configuration on tsunami propagation. Figure 17 shows the variation of normalized run-up heights (η/η₀) and normalized peak flow velocities (V/V₀) at simulated gauges (G3 and G4) with different tree arrangements, respectively, where η₀ and V₀ denote the run-up height and the peak flow velocity in the absence of vegetation. The results indicate that C. equisetifolia as a new vegetation layer at the back of existing P. odoratissimus (Collocation 3) can minimize the amplification of the peak flow velocity through the open gap, while generating a slightly larger peak flow velocity at the G3 station. The run-up height is higher than that of Collocation 1 (pure P. odoratissimus) and smaller than that of Collocation 2. This is because the peak flow velocity is reduced significantly due greater resistance of P. odoratissimus when it passes the P. odoratissimus vegetation layer, and in back zone of C. equisetifolia vegetation with smaller drag force, the velocity increment in open gap is smaller compared to that of other tree configurations.In Collocation 1 (pure P. odoratissimus), Figure 18 indicates that normalized run-up height and normalized peak flow velocity at G3 decreased monotonously with increased densities, while the normalized peak flow velocity of G4 first increased (vegetation density ≤ 0.15 tree/m²) and then presented adverse variation. Figure 19 indicates that the increased vegetation height and diameter at different growth stages of P. odoratissimus generate positively correlated effects on normalized peak flow velocity at G4, and significant attenuation effects on normalized run-up height and normalized peak flow velocity at G3. Generally, changing forest parameters (e.g. vegetation collocations, densities and growth stages) could be a more cost-effective way to mitigate tsunami damage.

4. Conclusions

Based on the shallow water equations of mass and momentum, an explicit depth-averaged 2D numerical model was developed to simulate wave attenuation and tsunami wave mitigation in a vegetated coastal region. The finite volume method was used to keep conservation of mass in this model. The proposed model was well verified for both breaking and non-breaking solitary wave propagations on a bare sloping beach, and long periodic wave propagation on a partially-vegetated sloping beach, as shown by comparison of experimental data. Overall, the present numerical model has been shown to accurately account for wave propagation and transformation in vegetated and non-vegetated sloping beaches.

In the case of full P. odoratissimus forest on an actual-scale beach, the water surface elevations at the front of the forest belt increased due to the effect of reflected waves caused by vegetation. Here, P. odoratissimus forest played a positive role in reducing maximum water surface elevations through the vegetated zone, and as such the maximum wave run-up height of the tsunami was effectively reduced. When a forest has an open gap, the negative effects can be reduced by choosing a beneficial arrangement of vegetation and open gaps. Sensitivity analysis was also carried out by changing tree collocation, tree densities, and growth stage, with results showing that they had important impacts in reducing the peak flow velocity of the tsunami wave at the open gap exit.

Overall, these study findings provide cost-effective natural strategies to improve the effectiveness of P. odoratissimus forests with an open gap by changing vegetation arrangements and vegetation parameters (collocations, tree densities, and growth stage), which are important in the design and establishment of vegetated bioshields against tsunami hazards.