{'fraction_non_alphanumeric': 0.026722207586933614, 'fraction_numerical': 0.012266201264488935, 'mean_word_length': 4.838973274370313, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Forest density is more effective than tree rigidity at reducing the onshore energy flux of tsunamis: Evidence from Large Eddy Simulations with Fluid-Structure Interactions

Communities around the world are increasingly interested in nature-based solutions to mitigation of coastal risks like coastal forests, but it remains unclear how much protective benefits vegetation provides, particularly in the limit of highly energetic flows after tsunami impact. The current study, using a three-dimensional incompressible computational fluid dynamics model with a fluid-structure interaction approach, aims to quantify how energy reflection and dissipation vary with different degrees of rigidity and vegetation density of a coastal forest. We represent tree trunks as cylinders and use the elastic modulus of hardwood trees such as pine or oak to characterize the rigidity of these cylinders. The numerical results show that energy reflection increases with rigidity only for a single cylinder. In the presence of multiple cylinders, the difference in energy reflection created by varying rigidity diminishes as the number of cylinders increases. Instead of rigidity, we find that the blockage area created by the presence of multiple tree trunks dominates energy reflection. As tree trunks are deformed by the hydrodynamic forces, they alter the flow field around them, causing turbulent kinetic energy generation in the wake region. As a consequence, trees dissipate flow energy, highlighting coastal forests reducing the onshore energy flux of tsunamis by means of both reflection and dissipation.


Introduction
The partial failure of a multi-billion dollar concrete sea wall upon the impact of the Tōhoku tsunami in 2011 invigorated the debate about possible nature-based approaches to mitigating tsunami risk (e.g., [1][2][3]). A few years before, damage assessments after the Indian Ocean Tsunami pointed to a reduction in wave impact near natural features [4][5][6][7][8][9][10][11][12], but the causal attribution of lessened damage to vegetation has remained contentious [13][14][15]. Adding to the challenge of clearly identifying the potential protective benefit of vegetation are field observations from other contexts. For example, the Japan Tsunami annihilated a coastal forest with an estimated 70,000 pine trees in Rikuzentakata, Iwate Prefecture, highlighting that there are no fail-safe approaches to countering a tsunami's destructive power. However, while the Rikuzentakata coastal forest was not able to withstand the 2011 event, it might have mitigated the impacts of the 1933 Showa Sanriku tsunami and the Chilean tsunami in 1960 [16,17].
One step towards adding clarity to the ongoing debate about the role of vegetation in tsunami-risk mitigation is to distinguish between different ecosystems and the differing effect they have on tsunami runup. Contrary to sea walls, vegetation cannot completely stop a tsunami, and its effectiveness depends both on the magnitude of the tsunami and on the structure and response of the vegetation. To leverage vegetation in tsunami-risk mitigation, we need to quantitatively understand the mechanism through which the vegetation interacts with the wave, the limits of this mechanism, and the degree to which it leads to a reduction in the energy of the tsunami [18]. The differences can appear quite subtle: For example, a detailed field assessment after the 1998 tsunami in Papua New Guinea found that Casuarina trees were able to withstand wave impact better than palm trees that were uprooted or snapped off at mid trunk level, transforming into projectiles and increasing rather than decreasing damage [19].
The lessons learned from both the Japan and the Papua New Guinea tsunami highlight that advancing our understanding of the protective benefit of coastal forests requires not only a quantification of the effect that the trees have on the tsunami but also an analysis of the effect that the tsunami has on the trees. One important difference between Casuarina trees and palm trees is the bendability of their trunks, with Casuarina trees having more flexible trunks. More bendability, however, is not necessarily better. In fact, field surveys in Sri Lanka and Thailand after the Indian Ocean Tsunami suggested that older and relatively sparse belts of Casuarina equisetifolia withstood the tsunami but did not appear to provide much protection [18], emphasizing the potential importance of forest density in addition to individual tree properties.
The goal of this study is to quantify how the energy dissipation and partial wave reflection of individual coastal trees vary with different degrees of bendability and the number of trees. We approach this problem numerically by solving the Navier-Stokes equations for turbulent flow with a free surface and fluid-structure interactions (FSI) using a multi-physics software package developed by some of the authors and described in [20,21]. The FSI module allows to capture the two-way coupling between the fluid flow and the mechanical bending of the trees due to wave forcing. While our work is motivated by prior field studies, we study the interactions between fluid flow and tree bendability in an idealized laboratory-like setting that is not intended to mimic any specific tree species or location. To enable transferability of our results, we consider a non-dimensional parameter range that is representative of tsunami flows and realistic hardwood properties [22].
Our approach of leveraging high-fidelity numerical methods complements existing approaches that adopted the depth-averaged shallow-water, Boussinesq or Serre-Green-Naghdi approximation in 1D or 2D (e.g. [23][24][25][26][27]). These models are powerful in capturing regional-scale inundation but are more limited in estimating forces on individual structures [28]. As demonstrated by comparing a 2D against a 3D model in [29], the solution of the 3D Navier-Stokes equations is required to quantify impact with the due precision. The limitation of not being able to capture 3D impact forces is less severe when the vegetation under consideration is relatively flat as compared to the water height, as might be the case with seagrass or a marshland, but is important in the context of bendable tree trunks experiencing tsunami impact.
Our work is not the first attempt to leverage Navier-Stokes-based simulations to improve our understanding of tsunami impacts on vegetation. Maza et al. [30] solved the Unsteady Reynolds Averaged Navier-Stokes (URANS) equations to study the interaction of a solitary wave with vegetation, where vegetation was modeled by rigid cylinders after the experiments by [31]. The same authors built a 2D linear model for vegetation motion and simulated the swaying dynamics of submerged vegetation [32]. An extension to 3D cylinders is presented in [33] where a rigid motion of the cylinders is allowed around a fixed point on the ground. Unlike what we show in this article, the above cited studies do not account for the cylinder bendability and elasticity along its length.
Navier-Stokes-based simulations of tsunami impacts come with their own challenges, such as the need for careful validation [34], sufficient resolution in boundary layers [35,36], and the potential overproduction of turbulence [37]. Fortunately, there is an increasing number of laboratory experiments that are essential for the validation of numerical simulations [12,[38][39][40][41]. Similarly, increasing computational power and data storage have solved some issues, particularly those related to insufficient resolution, making scientific computing for large scale tsunami simulations increasingly accepted [42,43]. This increased interest in high-performance computing for tsunami runup has led to several inter-comparison efforts published in recent years [44][45][46][47][48][49][50][51] that add valuable transparency to high-fidelity computational approaches.
While high-fidelity simulations provide a high level of detail regarding the flow field around submerged vegetation, post-processing is necessary to distill these details into a simple metric that is indicative of protective benefit. In much of the literature about the role of vegetation in coastal risk reduction, this metric is inundation distance, defined as the horizontal distance that a tsunami reached inland. In fact, one of the criticisms of the field studies that identified associations between vegetation and damage from satellite data [4,9] was that vegetation did not appear to have a mitigating effect on inundation distance [14]. However, inundation distance as a metric of protective benefit might be more valuable in the storm-surge context for which it was originally introduced [52,53] than for tsunamis. In contrast to storm surges, the main destructive potential of a tsunami is related to the energy transported not to the presence of water per se. We hence adopt reduction in energy flux as a metric of potential benefit, implying that there might be a significant potential benefit even if the inundation distance remains relatively unchanged, as seen in [23].
In this work, we implement a two-way coupling fluid-structure interaction model to investigate tsunami flow energy damping performance for different vegetation parameters. We use circular base cylinders as an idealized approximation for tree trunks following previous work on emergent, flexible plants in flow [54]. The remainder of this paper is organized as follows: § 2 describes the methodology followed by the model validation tests described in § 3. The numerical setup and the different scaling parameters are described in § 4. In § 5 we present the energy budget comparisons for different vegetation parameters, including velocity and turbulent kinetic energy distribution and wave profile evolution. In § 6, we discuss how our results could inform the ongoing debate on nature-based approaches to tsunami-risk reduction and end with a brief conclusion in § 7.

Numerical model description 2.1 Mathematical model
The interaction of a two phase fluid flow (water and air) with different sets of deformable structures is studied using a two-way-coupled fluid-structure interaction (FSI) model for multiphase turbulent flows. In order to track the moving interface between the solid body and the two phase fluid flow, referred to as the wet surface hereafter, the Navier-Stokes equations for large eddy simulation (LES) cast in arbitrary Lagrangian Eulerian (ALE) formulation are used [55]. The underlining idea of LES is the separation of scales in turbulent flows. This separation is usually achieved by a filtering operation that, for any vector or scalar quantity γ(x), is such that where the barred values are the filtered, resolved values which are separated from the un-resolved primed quantities at the sub-grid scale. The filtered Navier-Stokes equations of incompressible flows that account for the coupled deformation of solid structures are: where u i is the velocity vector of fluid i, p i is its pressure, ρ i is the fluid i's density, µ i is its dynamic viscosity, τ sgs i = (uu) i − (u u) i is the subgrid-scale stress tensor, g is the magnitude of the acceleration of gravity which points downwards in the direction of k = [0, 0, 1], and u m corresponds to the domain velocity which is obtained from the domain displacement given by the structure's deformation. The i subscript indicates the phase under consideration i = 1 for water and i = 2 for air. We model sub-grid scales using the approach of Vreman [56], described also in [57], with a model constant c vr = 0.1 used in all of the simulations.
At the free surface, the two fluids interact by means of a shear stress imposed via the two conditions: where t i represent the traction forces. We model the dynamics and deformation of the solid bodies representing the vegetation by the Euler equations as derived in [58]: where ρ s is the solid's density, d s represents the displacement field, P is the first Piola-Kirchhoff stress tensor, and b is the load distribution. At the wet surface Γ s , the kinematic and dynamic coupling conditions are imposed as: ∂t where d f is the displacement of the interface seen from the fluid perspective, t f and t s correspond to the traction forces exerted on the fluid-solid interface respectively, and u f is the velocity of the fluid. These conditions express continuity of the displacements, stresses, and velocities at the wet surface. We non-dimensionalize equation 2 using the length scale h avg , velocity scale u avg , and the stress scale ρ f u 2 avg [59] yielding: In our case, the length scale h avg represents the average dam bore height, and u avg corresponds to the average flow velocity of the bore before hitting the cylinder. This nondimensionalization results in non-dimensional parameters such as the density ratio (ρ r = ρ s /ρ f ), and non-dimensional modulus of elasticity E * = E/ρ f u 2 avg .

Numerical implementations
We discretize equations 1 and 2 in space via linear finite elements [21,60] and advanced in time via a backward differentiation formula of second order (BDF2) [61]. The fluid-structure interaction algorithm is detailed in [55]. For the solid mechanics problem, we discretize the Euler equations (2) using a standard Galerkin method for large deformations with a Newmark time integration scheme, details of the numerical strategy can be found in [62] and [63]. The FSI problem is solved with an iterative staggered multi-code approach. We use a strong Gauss-Seidel coupling algorithm with an Aitken under-relaxation factor. In each time step, we solve the fluid flow first to calculate the loads on the deformable structure required in the solid mechanics solver. We then solve the solid mechanics problem, and the new position of the body is sent back to the fluid mechanics solver. Convergence is reached when the residual of successive coupling iterations is less than the residual of coupling parameters. Details of the multi-code coupling strategy can be found in [55]. We use a conservative level set model to track the motion of the water-air interface. The level set equations are approximated by a stabilized finite element method using orthogonal subgrid scales for numerical stability as detailed in [64].

Boundary conditions
All boundaries (lateral, left, right and top walls) are specified as no-flux with free slip. We apply a no-slip boundary condition at the cylinder wall and at the bottom surface through the wall-modeled LES approach described in [65] with the following Reichardt's law of the wall: [66] u + = 1 κ ln(1 + κy + ) + 7.8(1 − e − y + 11 − y + 11 where u + = u uτ is the dimensionless velocity, y + = yuτ ν is the dimensionless wall distance, u τ = τw ρ is the friction velocity, τ w is the wall shear stress, ν is the kinematic viscosity and κ = 0.41 is the von Kárman constant.

Model validation
In the current study, we have selected two experimental analyses to validate the numerical model. We use Benchmark I mainly to validate our simple dam break computation without any obstacles. In contrast, Benchmark II is valuable to further assess our model in the presence of obstacles.

Benchmark I
In 1952, Martin and Moyce [67] executed a comprehensive set of measurements to characterize the dynamics of a collapsing liquid column on a rigid horizontal plane. We use the nondimensional data reported therein to validate our flow solver in the absence of obstacles. The choice of Benchmark I is particularly suited to validate a numerical model because, as observed by Martin and Moyce [67], the measurements were not influenced by either viscosity or surface tension and can hence be generalized across scales for both viscous and inviscid flow solvers. To compare the numerical results against the measurements, the following non-dimensional quantities are used: where η is the water depth and a is the initial base. The time evolution of the collapsing column at t * = 1.8553 and t * = 2.3191 is qualitatively compared in figure 1 against the experiments by Martin and Moyce whereas figure 2 shows the quantitative comparison of the front position versus the non-dimensional time.

Benchmark II
For Benchmark II, we compare the output of our code against the experimental data of [68].
The experiment consists of a flow generated by a collapsing column of water over a uniform, thin water layer. The interaction of the collapsing water column and the thin water layer is responsible for the formation of a bore propagating downstream. The domain is 16.6 m long,  0.6 m wide, and 0.45 m deep, and the water flume is separated in two parts. In the upstream region, the water is 0.25 m deep, and downstream the water depth is 0.02 m. The gate is situated at 5.9 m from the upstream wall, and the cylinder is placed 5.2 m from the dam. The diameter of the cylinder is 0.14 m. The schematic of the model is shown in figure 3. We have conducted our tests for different grid resolutions. Table 1 shows the details of grid sizes for fine and coarse resolution grid setup. The comparison of water depth evolution between the fine grid and coarse grid setup along the streamwise direction from the cylinder leading edge is given in figure 4. The location x c = 0 corresponds to the position of the leading edge of the cylinder, and x c is negative towards the upstream direction. At coarse grid setup, the maximum bore depth at cylinder leading edge reaches around 0.225 m at 3.6 s, whereas as reported by Arnason et al. [68], the maximum bore depth reached 0.27 m at 3.5 s, which indicates that our computation underestimates the water depth at this coarse grid resolution and is unable to capture the water level rise after the bore impact due to numerical dissipation.
In contrast, the fine grid resolution perfectly captures the maximum water level after the bore impact on the cylinder surface. Once the bore hits the cylinder at 3 s, the water level increases at the cylinder leading edge due to reflection, and the bore reaches a maximum of 0.268 m at 3.5 s. The details of the fine mesh grid set up and comparison of maximum water level at the cylinder upstream edge between the experimental and numerical case is reported in table 1. We also use the coarse and fine grid simulations to compare the depth-averaged flow velocity against the experimental data at the cylinder location (x = 11.1 m). As shown in figure 5, before the bore arrival at cylinder upstream edge, a negligible discrepancy of velocity exists between numerical analysis and laboratory experiment due to numerical oscillation  (t = 2 − 3.7 s), but once the bore reaches the cylinder upstream of the wall, the numerical result approaches to the average position of the velocity-time curve of the experiment. From the comparisons, we conclude that the current numerical method, including the turbulence model, wall function, grid size (0.0075 m), and the boundary conditions, correctly predict the water depth and flow velocity evolution and can be employed to simulate the tsunami bore and energy flux analysis for multiple cylinders configuration.

Numerical experiments
We carry out numerical experiments for different vegetation parameters such as elastic moduli, diameters, and stem spacing to calculate the energy balance. The details of numerical flume setup, non-dimensional parameters, and energy balance calculations are described in the following sections.

Size of the numerical flume
We modify our numerical setup used to validate our simulations against the experimental results by [68] to accommodate multiple cylinders in the domain along the spanwise direction. The modified testing flume is 20.6 m long, 0.8 m wide, and 0.45 m deep. The dam depth at the upstream reservoir is 0.18 m, and in the downstream region, the initial water depth is 0.02 m water depth. We increase the streamwise length of the domain to avoid wave reflection at the end wall and the domain width to minimize the blockage effect between the lateral walls and cylinders. The sketch of the numerical testing flume is shown in figure 6.

Scaled parameters of tsunami bore generation
The onshore height of a tsunami bore can vary from 2 m to 15 m based on the tsunami intensity scale [69][70][71][72][73]. In this study, we scale down the model to reduce the computational cost. Since the grid size of the numerical model needs to be sufficiently smaller to capture the small-scale eddies, the computational cost will be huge due to a large scale separation between the eddies unless the problem has been scaled down. Therefore, to simulate the tsunami phenomena, we scale down the computational model using the scale factor SF = 1 : 60. In this study, We choose a possible tsunami bore height of 4.5 m for the Indian Ocean tsunami [74] on a real scale, which corresponds to a bore height of 0.075 m in the scaled simulations. Following the method presented by [68], we then estimate the initial dam depth as 0.18 m based on the onshore tsunami depth in the scaled simulations. Our approach to scaling implies that the Froude number is approximately the same in a field setting as in our scaled simulations. We focus on maintaining scale similarity in the Froude number because it is an important parameter for onshore tsunami propagation and other freesurface flows, particularly where hydraulic jumps occur [75], [76]. A tsunami generally flows inland as a supercritical flow characterized by F r > 1. Different studies proposed different ranges of Froude numbers. For example, [77] suggested a range of 0.7 − 2 for the Froude number based on field surveys conducted for past tsunamis. Similarly, a study by [78] found Froude numbers between 1 and 1.5. Informed by these estimates, we use F r = 1.2 for all simulations, calculated based on the tsunami inflow condition before impact on the cylinders. We keep the inflow condition the same for all simulations, implying that F r is the same for all cases (see table 2).
The main drawback of ensuring scale-similarity in the Froude number is that the Reynolds in our simulations is lower than it would typically be in the field. Since the flow structure changes in the vicinity of the tree trunks due to dissipation, the Reynolds number (Re h ) is calculated for the scaled model using the hydraulic radius (R h ) as a characteristic length as employed by [68], where the hydraulic radius is defined as a ratio of flow area (A) and wetted perimeter (P ), which are defined as where h avg is the average incoming bore depth and B is the flume width. In this study, the Re h is 2.8 × 10 4 for scaled model, which represents the Re h value of 13 × 10 6 on a real scale.
In the scaled simulation, the initial bore depth is the same for all configurations. Since the hydraulic radius is calculated based on the flow area before the bore impact on vegetation, the value of Re h is the same for all configurations in scaled simulations. The bore parameters corresponding to the real scale and scaled model are presented in table 2.

Different parameters of multiple cylinder configurations
In these simulations, we have modeled trees as cylinders since emergent plants with rounded stems can be approximated as circular base cylinders [54]. Since we have scaled the tsunami flow parameters, the vegetation parameters must be scaled down as well. The main challenge is to scale the modulus of elasticity so that the deflection of cylinders will be correctly captured. Scaling the modulus of elasticity using the scale factor SF = 1 : 60 may not be correct because, apart from the modulus of elasticity, rigidity also depends on the moment of inertia, which can not be captured unless a proper scaling parameter is used. Hence, to properly scale the modulus of elasticity, the parameter EI, known as flexural rigidity, is scaled down. Here, the parameter I is the moment of inertia (I = πD 4 /4). A similar approach was implemented by [39], where rigidity was scaled using a non-dimensional vegetation flow parameter (V F P ) expressed as where H is the vegetation height, u avg is the flow velocity of the bore, B is the bed width, SP is the center to center distance between the cylinders, ρ r = ρ s /ρ f , ρ s is the vegetation density which is 500 kg/m 3 , ρ f is the fluid density (1000 kg/m 3 ) and E * = E/ρ f u 2 avg is the non-dimensional modulus of elasticity. Here, we choose hardwood trees such as pine as a reference point for our scaling. First, we calculate V F P based on the parameters on the real scale reported in table 2. The range of V F P used in this study is shown in table 3. Then, the scaled dimensional modulus of elasticity of the model can be calculated by the following  Since the rigidity (EI) changes with the diameter, we choose three sets of scaled moduli of elasticity and diameters for each case, totaling nine cases reported in table 4. The diameter of the cylinders is set in such a way that the effective channel width along the spanwise direction (Y direction) will be the same for all cases. Therefore, three cases with the same cylinder-perimeter are considered, one cylinder with a diameter of 0.14 m, two cylinders with a diameter of 0.07 m, and four cylinders with a diameter of 0.035 m. These configurations are chosen in such a way so that the rigidity (EI) is properly captured based on the vegetation parameters. Because once the diameter increases, rigidity also increases, and therefore to maintain the rigidity in the proper range, elastic modulus needs to be modified to model the cylinders properly. In this study, the range of diameters is chosen based on the rigidity of the scaled models used in the simulation. These configurations will provide more meaningful insight into what the best choice for tsunami mitigation should be, whether by increasing the elastic moduli or the number of cylinders with a smaller diameter. We also vary the gap between the cylinders to investigate how the spacing between the cylinders will affect the energy reflection and dissipation. Table 4 shows all the parameters of cylinders used in the scaled simulations.

Non-dimensional parameters
As stated in 4.3, V F P is used to scale down the modulus of elasticity. However, additional non-dimensional parameters need to be calculated to verify the scaled model to capture the tsunami bore propagation through vegetation. This set of non-dimensional parameters can either characterize the flow or vegetation configurations. In the review by Nepf [54], the nondimensional parameters such as the roughness density and solid volume fraction are introduced to measure the non-dimensional canopy density. As a starting point, and considering fully emergent vegetation, the non-dimensional measure of the vegetation density is the solid volume fraction [83].
where a = D/SP 2 is the frontal area per canopy volume for a given diameter D, and SP is the average distance between two consecutive stems. The parameters a and φ vary as a function of the tree diameter and spacing between the trunks. If the vegetation height is considered instead of the diameter, the non-dimensional canopy density is represented by the roughness density, which is defined as λ f = aH, where H is the canopy height. The range of both parameters  Table 4 shows the solid volume fraction and roughness density for each configuration. As reported by Mazda et al. [84], φ can reach up to 0.45 for dense mangrove forests. In the current configurations, since we are using pine forests as reference vegetation parameters, the value of φ should be lower than the solid volume fraction of dense mangrove forests. In our study, the solid volume fraction ranges from 0.2 to 0.35.
Since the tsunami becomes supercritical while propagating through the onshore vegetation, we calculate the Froude number before the bore impacts cylinders. As described by Arnason et al. [68], we calculate the surge Froude number, which is defined as F r s = u/ √ gh 0 where h 0 is the downstream water depth which is 0.02 m at t = 0 and u is the bore front celerity. The stationary Froude number at the cylinder upstream edge location is also calculated as v avg / gh avg where h avg is the average onshore bore depth and v avg is average flow velocity in scaled simulations (see table 2). Both parameters match well with the results reported by Arnason et al. [68], indicating that the flow becomes supercritical before hitting the cylinders. The Froude number, average water depth, and bore velocity for the individual cylinder configuration are compiled in table 5.

Energy Flux
The energy flux at a particular location gives information about how much energy is available for mixing and transport [85]. To quantify the protective benefits of tsunami mitigation, we calculate the energy flux at different streamwise locations as computed by [23,51]. The energy flux is defined as where q = u · u/2 is the kinetic energy per unit mass. Since, in this study, we aim to investigate the energy budget along the flow or streamwise direction, we integrate the energy flux at a fixed location x i along the water depth (along the Z axis) and spanwise direction (along the Y axis). We then normalize the depth and spanwise integrated energy flux by the width of the domain so that the energy flux will only vary along the horizontal, streamwise direction.
The cumulative time-integrated horizontal energy flux is defined as where F t is cumulative time-integrated total energy flux. We take two transects along the streamwise direction to calculate the reflected horizontal energy flux. Transect T-I is taken in the upstream region (x = 8 m) and another transect T-II is near to the cylinder leading edge (x = 11 m) as shown in figure 7. If there is any reflection, the energy flux will be drastically reduced at T-I due to the reflected wave arrival. We integrate the horizontal energy flux over the flow period (see equation 13) at T-I and T-II to identify reflection. Finally, we calculate the reflected cumulative horizontal energy flux by taking the difference between the cumulative horizontal energy flux at T-I and T-II.
where F k R and F k R are the reflected cumulative horizontal kinetic and potential energy flux, respectively. Finally, we obtain the reflected energy flux coefficient by normalizing kinetic and potential energy flux by the time-integrated total energy flux at T-I, which is calculated as: We calculate the cumulative energy flux at T-III, which is located behind the cylinder trailing edge, and at T-IV, located at 2 m from the cylinder trailing edge (see figure 7). Then, the dissipated cumulative horizontal energy flux is calculated by subtracting the cumulative energy flux at the T-III and T-IV, as shown below.
where F k D and F k D are the cumulative dissipated kinetic and potential energy flux, respectively. Similarly, the energy flux dissipation coefficient is defined as: Since the bore generated in the wet bed reaches the transects (T-I, T-II, T-III, T-IV) Table 6 summarizes the energy flux parameters calculated in this study.

Rigidity enhances energy reflection for a single cylinder
Energy reflection is one of the primary means of attenuating tsunami waves [23]. When the bore hits the cylinder, part of its energy is reflected. To identify the reflection, we show the free surface profile at different flow times for one cylinder in figure 8. At t = 4 s the bore front reaches the cylinder upstream face with a bore depth of approximately 0.075 m. As the bore continues to propagate, the bore depth at the cylinder's upstream face continues to increase due to flow blockage and reaches a maximum of 0.14 m before reflection begins. At the downstream side, the split flow joins and begins to form a turbulent wake, whereas the reflected wave flattens out towards the sides of the flume and continues to propagate back upstream (see figure 8b). To quantify energy reflection, we calculate the energy flux reflection coefficients as described in section 4.5. To study how reflection varies with the vegetation properties, we compute the reflected energy flux coefficients for different elastic moduli. Figure 9 shows the comparison of   As stated in section 4.5, the reflected energy flux coefficient is calculated based on the cumulative time-integrated reflected energy flux, which means that the slope of the reflected energy flux coefficient shown in figures 9 indicates the rate of change of reflected energy over the flow period. During the initial bore impact and propagation (t = 4.4 − 7 s), the slope of the curve of reflected kinetic energy flux coefficient becomes steeper for the cylinder with the highest modulus of elasticity than with the lowest modulus of elasticity, indicating that the amount of reflected energy increases as rigidity increases. As the bore propagates further downstream, the gradient of reflected kinetic energy flux coefficient becomes constant at higher rigidity. However, at the lowest modulus of elasticity, the gradient still increases (see figure 9) for t = 8−10 s. Since the rigid cylinder deforms less after the bore impact but also rebounds less with a negligible deflection after passage of the bore front, the reflected energy is higher during the initial bore impact, and becomes constant quickly. On the other hand, flexible cylinder deforms more by the hydrodynamic forces exerted on them and hence reflects less energy. However, the flexible cylinder also rebounds more compared to the rigid one, resulting in a continuous increase of reflected energy during the later phase of flow propagation (t = 8 − 10 s).
The reflected potential energy flux coefficient is consistent with the evolution of the reflected kinetic flux coefficient. In figure 9 we find the gradient of reflected potential energy flux coefficient is steeper at the highest rigidity, but during the later phase of flow, the gradient declines for all cases, although the slope declines more rapidly as the rigidity increases. Since rigidity triggers reflection, the bore depth increases at the upstream face of the cylinder. This finding is supported by the average water surface evolution at the cylinder upstream edge location at different elastic moduli shown in figure 10. Since the potential flux depends on the water depth, rigidity-induced water elevation increases the potential flux.

Reflective effect of rigidity is lost for multiple cylinders
Rigidity is not the only parameter affecting onshore energy flux; the number of cylinders and the cylinder spacing also alter the flow field sensitively and hence the energy flux. Interestingly, our simulations show that the dependency of reflected energy on the rigidity is most pronounced for the single cylinder shown in figure 9. The equivalent results for multiple cylinders are summarized in figures 11a and 11b. In table 7, we compare the percentage of maximum In a two-cylinder arrangement, this difference decreases to 2%. In a four-cylinder arrangement, this difference further decreases to 1%, suggesting that the effect of rigidity on the reflected kinetic energy flux becomes negligible as the number of cylinders increases. In a multiple-cylinder configuration, the blockage area along the span-wise direction is split into smaller areas. The flow can propagate through the gap between the cylinders, leading to a different flow field with multiple wakes. Figures 12 and 13 show two snapshots of the flow around two cylinders and four cylinders at the time of maximum deflection and after passage of the bore front, respectively. At t = 4.2 s after bore impact, the increase in water depth at the cylinders' upstream faces is clearly visible for two cylinders. The interactions between the four cylinders are apparent, particularly in figure 13, where the split flow around each one of the cylinders joins into one wake that is reminiscent of the wake behind a single cylinder in figure 8 (b).
For four cylinders, these flow-field interactions become even more pronounced, as shown in figure 13. In this configuration, the cylinders no longer deform in unison, because of the dynamic feedback between deflection and flow. In figure 13a only the second cylinder from the upper side deflects primarily in the flow direction while the cylinder closest to the lower side of the domain bends partially outwards to accommodate flow in the gap between the cylinders. As a consequence. the reflected wave front in the four-cylinder configuration does not increase in depth, but instead flattens out along a spanwise direction. However, the total blockage area of two cylinders is twice of four cylinders, which explains why increasing the number of split regions between the cylinders reduces the energy reflection. The dynamic feedback between cylinder deflection and flow explains why the relationship between gap and energy flux is not a simple one. Our simulations demonstrate that the gap between the cylinders alters flow restriction and modifies the energy flux reflection, but the difference between different configurations is relatively small. The blockage area of cylinders are the same at SP/D = 1.5 and SP/D = 2; as the gaps are widened, the cylinders are shifting towards the sidewall, which explains why increasing the gap does not change the energy reflection much. As shown in table 7, the difference between SP/D = 1.5 and SP/D = 2 in energy flux reflection is 0.3% for two-cylinder and 0.2% for four-cylinder configurations and hence negligible for practical purposes. Since the blockage area of a two-cylinder configuration is two times of that in a four-cylinder configuration, the cylinders with a higher blockage area cause slightly higher reflection in closely spaced cylinders.
The percentage of reflected kinetic energy flux lies between 25 − 30% for all configurations, whereas the amount of maximum potential energy flux reflection is 4 − 5 times smaller than the kinetic energy flux reflection. The reflected potential flux shows the same pattern as kinetic flux, where the importance of rigidity is lost when the number of cylinders increases. Since potential flux depends on the water depth at the cylinder upstream edge (transect II, see figure 7), the comparison of water depth evolution for different elastic moduli is shown for multi-cylinder arrangement in figure 14 to verify out findings. Since the difference in energy flux reflection at SP/D = 1.5 and at SP/D = 2 is minimal (see table 7), the water surface evolution for different flexibility is shown only at SP/D = 1.5. Also, the difference in wave height for different rigidities at the cylinder upstream edge becomes minimal in a two and four-cylinder arrangement, indicating that the flexibility-induced wave reflection is negligible in multi-cylinder configurations.

Flow blockage area enhances the wave reflection
In this study, we define three different cylinder configurations, each cylinder configuration has a different flow blockage area. The flow blockage area in a one-cylinder is two times higher than in a two-cylinder and four times higher than in a four-cylinder configuration. As shown in table 7, the one-cylinder with the highest flow blockage causes more flow restriction than multiple cylinders, resulting in the maximum amount of energy reflection. To compare the wave reflection for all cylinder arrangements, the free-surface evolution for different cylinder arrangements at transect T-II (x = 11 m) is shown in figure 15. Without the cylinders, the bore maintains a constant depth over the flow period. The bore depth increases when cylinders are present, indicating the water level rise due to flow blockage, which causes the reflection of the wave. The reflected bore depth is highest in one cylinder with a higher flow blockage area and higher in closely spaced cylinders. To compare the wave reflection for different cylinder configurations, we measure the maximum reflected bore depth (H R ) by taking the maximum    value of wave height over the flow period. We define the incident bore depth in the absence of cylinders as the maximum bore depth (see figure 15). Finally, we obtain the reflection index by taking the ratio of H R and H B and comparing it for different cylinder configurations. A similar approach was implemented by [86] to calculate the reflection index of tsunami bore to study the influence of different onshore seawalls on tsunami wave attenuation. From the reflection index shown in figure 15, we conclude that the reflection of the wave follows the same pattern as the energy flux reflection variation in different cylinder arrangements. As expected, the reflection index in one cylinder with a larger diameter is highest among the other configurations, whereas the four-cylinder with a larger gap (SP/D = 2) has the minimum reflection index. The percentage of reflection for different cylinder configurations is also shown in table 8. The smaller gap between the cylinders is beneficial in terms of the reflection because cylinders with a smaller gap produce more flow blockage. Nevertheless, when the number of cylinders is increased to four, the effect of closely spaced cylinders in flow blockage becomes negligible. The percentage of wave reflection for the four-cylinder with a smaller gap (SP/D = 1.5) is only 1% higher than the four-cylinder at SP/D = 2.

Rigidity reduces flow damping and turbulent kinetic energy at downstream region
Understanding the flow characteristics in the downstream region is important because it gives us an idea of how the flow field changes in multi-cylinder configuration and at different rigidities.
In the wake region, the bendability not only controls the motion of the cylinder body but also flow structures around the. Since not all cylinders move uniformly, the dynamic response between the flow and the cylinder bodies makes the flow structure more complex, which affects the flow velocity and hence can alter dissipation. To understand the effect of flexibility in the wake region, we show the longitudinal average velocity profile along the streamwise direction in figure 16. The velocity is averaged along the depth and spanwise direction. The field averaged velocity is then further averaged over the flow period after the bore has reached the downstream region. Here, the flow period is 4 s which is the time interval between the bore arrival at the downstream region and the end wall (see figures 8, 12 and 13). The velocity is normalized by the approach velocity (U ∞ ), which is the velocity of the bore before hitting the cylinders to compare how the longitudinal velocity changes along the wake region with respect to the approach velocity. As shown in figure 16, the averaged velocity is lower than the approach velocity at the upstream edge of the cylinders. Once the bore flows past the cylinders, the flow decelerates, and the flow energy is reduced due to reflection. At the flow separation region, the normalized velocity steadily increases from 0.42 and maintains a constant value at a further downstream region. As expected, flexibility induces more velocity reduction in the wake region. The normalized velocity is reduced by 15% at the streamwise location x = 11.4 − 12 m. At the farthest downstream region (x = 11.8 − 12 m), the normalized velocity is reduced by around 20% at modulus of elasticity E = 5 × 10 5 Pa in four-cylinder arrangement indicating that flexibility alters the flow structure and the change of the flow structure at the downstream region increases as the number of cylinder increases. Figure 16 shows the longitudinal velocity profiles at SP/D = 2 only because we have not observed any significant changes in normalized velocity as the gap between the cylinder changes. The longitudinal profiles of instantaneous turbulent kinetic energy (TKE) in figures 17 and 18 show the effect of flexibility and gap between cylinders on the reconfiguration of turbulent flow structures in the wake region of cylinders. The turbulent kinetic energy per unit mass is calculated by taking a spanwise and depth average over the velocity fluctuation fields as: where u = u − u , v = v − v and w = w − w are the velocity fluctuations and represents the spanwise average operator. Finally, the turbulent kinetic energy is normalized by k ∞ measured at x = 11.1 m which is the location of the upstream cylinder edge. For all the investigated flows, the normalized instantaneous field-averaged turbulent kinetic energy profile (k/k ∞ = 1) is extracted at t = 6 s, representing the scenario when the dam front wave has reached the downstream region (see figures 13 and 12).  18 show two peaks of turbulent kinetic energy: the first peak is located behind the cylinders and is a consequence of flow separation. The second peak occurs further downstream, where the split flow joins, and a hydraulic jump occurs. For both the twocylinder and the four-cylinder arrangements, the first peak of normalized TKE becomes higher as flexibility increases at the trailing edge of cylinders. We detect the highest first peak at SP/D = 2, where the normalized TKE increases by 30% at the lowest modulus of elasticity. After the first peak, a region of reduced TKE extends up to x = 11.8 m for all scenarios. Between x = 11.8 m and x = 12 m, a second peak forms, initiating a von Kaŕmań vortex street. The highest second peak arises at SP/D = 2 for both two and four-cylinder arrangements at the lowest modulus of elasticity. For four cylinders, the difference between TKE at E = 5 × 10 5 Pa and E = 1 × 10 6 Pa is higher at SP/D = 2 than at SP/D = 1.5. A similar pattern emerges in the two-cylinder arrangement. The reconfiguration of the cylinders alters the flow field and enhances the turbulence generation in the downstream region. Also, the flexibility alters the turbulent structures in the wake of cylinders, particularly at SP/D = 2.

Rigidity reduces kinetic energy dissipation for all configurations
Energy reflection is not the only physical mechanism that contributes to reducing the energy carried onshore by a tsunami. Energy dissipation can help reduce the energy flux created by the tsunami and is hence an important parameter to assess the protective benefit afforded by different vegetation configurations in tsunami mitigation. Figure 19 is the equivalent of figure  9, showing the dissipated rather than the reflected energy flux coefficient for a single cylinder. It demonstrates that increasing rigidity has opposing effects on the dissipation of kinetic and potential energy fluxes: More rigidity dissipates less kinetic energy but more potential energy. Once the bore hits the cylinders, the flexible cylinder deforms more than the rigid cylinder, dissipating more kinetic energy. After the impact of the bore front, the energy of the incoming flow is not sufficient to cause further deformation of the cylinder. Therefore, the dissipation rate of kinetic energy becomes almost constant at the later phase of the flow period (t = 8 − 10 s) in figure 19. The dissipated potential energy flux is around 80% of the dissipated kinetic energy flux coefficient in the one-cylinder configuration. The effect of rigidity in the potential flux Figure 19: Energy dissipation flux coefficient as a function of time of one cylinder for different elastic moduli. φ k d dissipated kinetic energy flux coefficient. φ p d : dissipated potential energy flux coefficient dissipation is exactly the opposite of what we observe in kinetic energy dissipation. The highest potential flux is observed in a less deforming cylinder. The more pronounced deformation of the flexible as compared to the rigid cylinder also reduces build-up of the free surface around the cylinder and hence dissipates less potential energy flux. As a consequence, the effect of rigidity on the total dissipation, including both the kinetic and the potential component, is small.
We observe the same pattern of energy dissipation in multi-cylinder configurations. Figure  20 shows that rigidity continues to decrease the dissipation of kinetic energy but increases the potential flux dissipation even in the presence of multiple cylinders. The tendency for rigidity to increase the dissipation of potential energy flux remains notable for two-cylinders in figures 20a and 20b, but disappears in the four-cylinder configuration (20c, 20d). Since we do not observe any noticeable build-up of the wave surface around the cylinder, the variability of wave depth along the downstream region is small in multi-cylinder as compared to onecylinder. The dissipated potential energy flux becomes 50% of dissipated kinetic flux in the two-cylinder configuration. In the four-cylinder configuration, the dissipated potential flux is further reduced to 30% of kinetic flux. Also, with increasing rigidity, the dissipated potential flux increases. Since flexibility alters the wave reflection, the difference in wave height between the upstream and downstream regions becomes smaller. However, as the number of cylinders increases, flexibility has less impact on dissipated potential energy flux. In the two-cylinder configuration, the dissipated potential flux decreases by 5% as flexibility increases. In a fourcylinder configuration, the changes of potential flux are less than 1% since the changes in the depth of the transmitted wave are negligible as flexibility changes.
We report the percentage of maximum dissipated kinetic energy flux for different cylinder parameters in table 9. For one cylinder, the percentage of the maximum dissipated kinetic flux coefficient is around 24% at an elastic modulus of 1 × 10 5 Pa, and the dissipated kinetic energy increases by 4% as the structure becomes more flexible. In the two-cylinder and four-cylinder configurations, the percentage of dissipated kinetic energy flux coefficient for the lowest modulus of elasticity steadily goes up to 33%. For both scenarios, higher flexibility induces more energy dissipation. Moreover, the difference in kinetic energy flux dissipation increases as the number of cylinders increases. In a two-cylinder configuration, the difference is 4%, whereas, in a fourcylinder, the difference becomes 5%. Our finding is also supported by the turbulent kinetic energy shown in figures 18 and 18, which prove that flexibility-induced dissipation enhances the Top Row: dissipated kinetic and potential energy flux coefficients for two-cylinder configuration. Bottom Row: dissipated kinetic and potential energy flux coefficients for four-cylinder arrangement turbulent kinetic energy due to dissipation. The comparison shows that the energy dissipation increases with an increasing number of cylinders, suggesting that the structure with a larger diameter contributes to energy reduction by reflection primarily. The contribution of potential flux dissipation is smaller in multi-cylinder configuration than in one-cylinder. The total energy flux dissipation is maximum in the four-cylinder configuration. However, the difference in total dissipated energy flux is not significant as the number of cylinders increases. Apart from flexibility, the gap between the cylinders also modifies the energy dissipation. However, the gap between the cylinders has a minimal effect on energy dissipation, the same phenomenon we notice in energy reflection. The difference in kinetic energy dissipation is 2% − 3%for all configurations. At the highest modulus of elasticity, where cylinders are almost rigid, the maximum energy dissipation obtained is approximately 2.2% higher in the four-cylinder at SP/D = 2 than at SP/D = 1.5. Whereas, In the two-cylinder arrangement, the difference is small (0.6%) for the rigid body, which indicates that increasing the cylinders alters the energy dissipation.

Discussion
The hardening of the shoreline by way of concrete sea walls comes at a staggering price, not only in terms of the construction costs, but also in terms of long-term negative impact on coastal ecosystems [87][88][89] and shoreline stability [90,91]. Less obvious than the drawbacks of lining the coast with concrete, often meter-high sea walls is the alternative, particularly along high-impact coasts like Japan or Indonesia, where tsunamis are frequent and destructive. The idea that vegetation may act as a bioshield against flooding was first proposed in the context of marshlands attenuating storm surges in southern Louisiana (e.g., [52,53]), but the wrath of a tsunami differs fundamentally from the comparatively mild tidal or seasonal flows for which vegetation is thought to be an effective bioshield [92]. An important aspect of tsunami runup is its profound dependence on both off-shore bathymetry and on-shore topography [93][94][95][96][97][98]. While off-shore bathymetry is difficult to modify, on-shore topography could be altered in a way that reduces tsunami impacts, essentially combining nature-based elements and traditional engineering elements into a hybrid approach to mitigating tsunami risk. An example of this approach is a vegetated coastal mitigation park: A landscape unit on the shoreline built specifically to protect communities or critical infrastructure and provide vertical evacuation space while strategically integrating trees or other vegetated elements to enhance protective benefit.
Prior work has characterized the protective benefit of the onshore topography created by the hillscape of coastal mitigation parks [23], but did not explicitly consider the role of vegetation. As communities across the world increasingly consider or adopt coastal mitigation parks, it remains an open question whether the presence of forests of large trees such as the maritime pines of Japan could be used to enhance the protective benefits of coastal mitigation parks and, if so, whether some species of vegetation might be more effective than others. The latter question has immediate relevance for ongoing international efforts to mitigate tsunami-risks. For example, field assessments of damage associated with the 2004 Indian Ocean earthquake and tsunami have posited that Casuarina forests mitigated tsunami impacts more effectively than native ecosystems like sand dunes [99]. As a consequence, natural sand-dune ecosystems along the Coromandel Coast, India, were destroyed and replaced with Casuarina Equisetifolia, an exotic timber with potentially adverse ecological effects [100].
Several field assessments have argued that different species vary significantly in the protective benefits they provide, as reviewed in [1]. However, many of these assessments are based on observed correlations between tsunami damage and the presence of different ecosystems [4][5][6][7][8][9][10][11][12]18] and stop short of identifying a possible causal relationship between vegetation and tsunami damage. A key challenge in comparing the effect of different ecosystems is that they thrive in different environments. For example, after the Indian Ocean tsunami Chatenoux and Peduzzi found that tsunami damage was reduced behind seagrass beds but increased behind coral reefs [9]. As pointed out by the authors, an important confounding factor is the difference in the morphological settings where these ecosystems are located with seagrass meadows growing primarily on shallow, gently-sloping continental shelves and barrier reefs lining a break in the continental slope.
Here, we intentionally remove the complexities associated with the geomorphological setting and consider a highly simplified flat topography to isolate the effect of vegetation properties. Since we are interested primarily in the role of onshore vegetation in tsunami-risk reduction, we focus specifically on coastal forests and the rigidity of tree trunks. In 5.1, we show that rigidity enhances energy reflection for a single cylinder, but the case of a single cylinder is primarily of abstract value since tree trunks rarely grow in isolation. In the practically more relevant case of multiple cylinders in 5.2, the difference created by variable rigidity all but disappears. Moreover, in 5.3, we observe that the gap between the cylinders has less effect on energy reflection as the number of cylinders increases. The additional complexity introduced into the flow field by the presence of multiple cylinders dominates over the rigidity effect, even in the idealized setup of our numerical experiments. There would be significantly more noise in an actual field setting, further reducing the difference in response.
Taken together, our simulations from sections 5.1 and 5.2 suggest that the rigidity of tree trunks is less relevant than the arrangement and density of the tree trunks within the coastal forest. An additional factor that reduces the role of rigidity is that rigidity reduces the dissipation of kinetic energy flux afforded by the presence of trees. At least for a single cylinder, rigidity hence affects the reflection and the dissipation of the onshore energy flux in opposite ways. Compared to the effect of rigidity on the reflection of kinetic energy flux, the reduction in dissipation of kinetic energy appears less sensitive to the presence of other tree trunks as demonstrated in section 5.5. However, the effect itself is on the order of a few percent and may hence not generalize to a field setting where many other factors confound flow behavior.
An interesting aspect of our findings in sections 5.1 -5.4 is that tree trunks reduce the onshore energy flux to a comparable degree through both reflection and dissipation. In contrast, the protective benefit of coastal mitigation parks lies primarily in the reflection of incoming wave energy [23]. This finding suggests that the protective benefit of tree trunks and hillscapes could complement the protective benefit of coastal mitigation parks by augmenting reflection and adding dissipation. An important caveat to this interpretation is that we only consider reversible, elastic deformation of the tree trunks and do not include tree breakage or uprooting, which tend to be the main cause for trees increasing damage in the field (e.g., [101]). The breaking patterns of trees do tend to vary among different species [102] and also depend on trunk diameter [103]. Particularly the dependence on diameter highlights the value of working with existing, native vegetation, and potentially integrating it into mitigation park, rather than planting a mono-culture of trees of a particular species from scratch.
We emphasize that our work adopts a highly idealized view and does not capture sedimentation and erosion that could significantly affect the stability of tree trunks. Neither do we consider the role of the canopy or other vegetation that would undoubtedly grow around trees. By reducing trees to cylinders with a constant elastic modulus, we neglect possible speciesspecific differences in the breaking pattern of the trunk, the branching structure emerging from the trunk or the root structure. Our work is hence merely a small step towards improving our understanding of the intricate interactions between onshore vegetation and tsunami impacts.
Further work would be needed to assess the protective benefit of vegetation with an inherently three-dimensional structure like mangroves that can not be meaningfully reduced to a single trunk or stem.

Conclusions
In this study, we quantify how the reflection and dissipation of onshore tsunami flow depends on the rigidity and arrangement of trees in a coastal forest. We use a three-dimensional LES turbulent model to study the effect of varying rigidity, vegetation gap, and the number of trees on tsunami energy attenuation. To represent the dynamic interplay between the tsunami bore and the trees, we employ a two-way FSI model to capture the coupling between fluid and solid body interaction. We compute energy flux at different streamwise locations to quantify the energy reflection and dissipation for three sets of moduli of elasticity and different cylinder configurations. Our simulations suggest that the rigidity of tree trunks only increase the reflection of kinetic energy flux in the somewhat abstract case of a single cylinder. Even small deflections of the tree trunk alter the flow structures downstream, impacting flow velocity and enhancing turbulent kinetic energy. This variability in velocity distribution and turbulent kinetic energy is more pronounced when the number of cylinders increases, highlighting that rigidity is less important than the number and arrangement of the trees. An important caveat to this conclusion is that we did not capture shear-driven erosion, sediment transport or tree breakage, all of which could limit the ability of the tree trunks to withstand tsunami impact.