Multiscale evaluation method of the drag effect on shallow water flow through coastal forests based on 3D numerical
 simulations

This study presents a method for determining the drag parameter in the 2D shallow water (SW) equation for flows through a coastal forest by conducting a series of 3D numerical simulations (3D NSs). Following the theory of multiscale modeling, an evaluation method procedure is proposed. We first prepare a local test domain that contains a sufficient number of trees to constitute part of a coastal forest. Then, 3D NSs are conducted in this test domain with various inflow conditions. Based on the corresponding results, the momentum losses over the test domain are converted into the drag parameter of the global SW equation. A response surface of the drag parameter is constructed as a function of the flow conditions. The stabilized finite element method is employed for both the local and the global NSs, and the phase‐field method is utilized to represent 3D free surfaces. Comparisons between the 2D SW calculation results and the 3D NS results are also performed to verify the validity of the proposed method.


INTRODUCTION
Coastal forests have received considerable attention in recent years as an ecodisaster risk reduction (Eco-DRR) management system for water related disasters. Since Shuto 1 published his pioneering work, many researchers have investigated the effectiveness of employing forests to mitigate the effects of coastal tsunamis. 2,3 Some observational investigations on the 2004 Indian Ocean tsunami reported that areas with mangroves trees were less damaged than areas without them. [4][5][6][7] The presence of black pine trees planted along the northeastern coast of Japan is also assumed to contribute to attenuating the power of the 2011 tsunami. 8 Comprehensive and in-advance assessments with numerical simulations (NSs) are a practical strategy for future tsunamis.
In previous shallow water (SW) tsunami/flood simulations, 9,10 the presence of vegetation has generally been taken into account by increasing the Manning roughness. According to the comprehensive review in Bricker et al., 11 many previous studies typically set a relatively high roughness parameter in forest regions based on early observational works; see, for example, Kotani et al. 12 and Te Chow. 13 In this context, an equivalent roughness model was developed by Petryk and Bosmajian 14 as an advanced estimate, and this model has been widely employed in NSs of tsunamis 15,16 and floods. 17 However, the original concept of roughness modeling is based on the idea of bottom frictional stress, which means that the resistance stress should decrease with increasing depth, as confirmed in an observational study. 18 The resistance of an emerged structure is assumed to increase as the water level increases. Thus, Bricker et al. 11 noted that such a discrepancy could be a source of uncertainty in the whole simulation. To resolve this issue, an alternative idea was proposed by Tanaka et al. 19,20 for coastal trees. This formulation addresses the drag effect of trees in terms of the flow velocity, depth, and drag coefficient. Since the resistance force increases as the depth increases, current studies that include not only water flow analysis but also wind flow analysis 21 are likely to employ this type of drag force modeling for trees. In addition, some research has been performed on the resistance characteristics of vegetation by either numerical analysis, laboratory experiments, or theoretical derivations. [22][23][24][25] With the aim of effectively evaluating the drag effect in the aforementioned approaches, multiscale modeling has attracted attention. In this type of modeling, the macroscopic damping effect of trees on the global flow is associated with the cell-averaged values of the local flow characteristics. [26][27][28] As emphasized in Liu et al., 28 numerical analyses based on multiscale modeling have a computational efficiency advantage over 3D simulations, and they provide deeper physical insights than macroscopic empirical approaches. Based on the principles of this approach, predictions of global flow through coastal forests can be efficiently conducted without ignoring the local flow characteristics around individual trees.
In this study, following the theories of multiscale modeling, we propose a method for determining the drag parameter in 2D SW flows within a coastal forest at the global level by conducting a series of 3D NSs of water flowing through a group of trees at the local level. The surrogate modeling employed for the global drag parameter evaluation, which is a well-known concept in machine learning or statistical analysis, is also one of the original elements of the proposed framework. The proposed method is divided into four processes. First, we prepare a local test domain (LTD) that contains a sufficient number of trees to constitute a coastal forest. Second, a series of 3D flow simulations are conducted in the test domain with various inflow conditions. From these simulations, representative responses can be obtained, which depend on the local flow conditions, flow velocities, or depths. Third, on the basis of each of these local responses, the pressure loss over the local domain is converted into the drag parameter of the global SW equation. Finally, a response surface of the drag parameter is constructed as a function of both the inflow velocities and the water depths. By virtue of the 3D NSs conducted at the local level, the complex flow characteristics associated with tree canopies are reflected in the drag parameter for the global SW flow. The resulting response surface is a function that acts as a surrogate for a number of time-consuming 3D NSs at the local level.
To verify the validity of the proposed method, we conduct 2D SW simulations, in which rigid tree models are located in an open channel, by using the identified function of the drag parameter. Comparisons between the 2D SW calculation results and the 3D NS results are also performed. The presented verification analysis exemplifying the concrete procedure simultaneously demonstrates that the proposed method enables a reasonable evaluation of the disaster mitigation effects of coastal forests with affordable computational costs without conducting 3D full-scale NSs at the global level. Although the presented framework is intended to contribute to improving the quality of the simulations of 2D SW flow through trees, it neither focuses on the characteristics of specific plant species (black pines, mangroves, etc.) nor investigates the possibility of incorporating arbitrary previous tsunami/flood models. For both the local and global NSs, we employ the stabilized finite element method (FEM) and incorporate the phase-field method to represent free surfaces in 3D NSs.

Summary and procedure
Our evaluation of macroscopic tree drag effects is based on the idea of spatial-scale separation, which distinguishes local and global flows. The characteristics of the local flow reflect those of the global 2D flow by a spatial averaging procedure along the lines of the multiscale theories. The proposed evaluation method consists of the following four steps and is illustrated in Figure 1: (i) We prepare a "LTD" that contains a sufficient number of trees in a coastal forest, as shown in Figure 1A. (ii) With various flow conditions at the local level, as illustrated in Figure 1B, a series of 3D NSs are conducted in a rectangular open channel, where the LTD is equipped. (iii) On the basis of each of the numerical test results, the macroscopic pressure loss relevant to the overall resistance of the trees is converted into the drag parameter for the global flow, and the flow velocity and depth are calculated from the volume-averaging operation over the LTD; see Figure 1C. (iv) The evaluation results are transformed into a response surface representing the drag parameter in the global flow as a function of the input parameters given as local flow conditions for the numerical tests.
The volume-averaging processes is performed over the LTD, which is regarded as a representative elementary volume (REV) 29 within the framework of porous media theory. Additionally, in this study, the 3D NSs in (ii) are referred to as "numerical flow tests" in the context of multiscale computational homogenization. 30 Note that the numerical flow tests are performed for the entire domain of the open channel, which contains the LTD, as illustrated in Figure 1B,C.
In the following section, each process is explained in detail along with the underlying theories.

Setup of the LTD
The first step of the procedure is to set up the LTD from the whole structure of a coastal forest that contains a sufficient number of trees to be representative of the homogeneous feature. The extracted LTD is a spatially reduced subdomain of the overall forests and is equipped in a rectangular channel, as shown in Figure 1A. This subdomain is used for the 3D NSs to numerically "measure" the trees' resistance to flow in the next step.
In some previous studies with purposes similar to ours, the length of the LTD in the overall flow direction was reduced to accommodate only a few trees. For example, Wang et al. 27 employed a periodic LTD, called a unit cell, in which cylinders are vertically located to develop a surface water model. Lee and Yang 31 also investigated flows in a 2D unit cell with emerged cylinders. However, Nepf et al. 32 reported that the drag induced by downstream trees is less than that induced by upstream trees. Their insights imply that the periodic unit cell in the streamwise direction is not appropriate to evaluate the total decay observed in the overall forests. In adopting this knowledge, we set the length of the LTD in the overall flow direction to L, which is the same size as that of the whole domain. Furthermore, the width of the LTD, Δw, can be set to be considerably smaller than the width of the overall forest, w, to contain only a few trees since the transverse direction has little effect on damping. Maza et al. 22 reported that the distribution of trees, which is quite random in a native forest, is not as influential as the density. Therefore, trees are assumed to be homogeneously arrayed in the LTD to have an equivalent population per unit area to the original forest. Specifically, the arrangement can be periodic, as often assumed in various approaches for multiscale modeling.

Numerical flow test in the LTD
With the LTD introduced in the previous section, 3D NSs are conducted as numerical flow tests to characterize the drag effect of a coastal forest. As mentioned in the previous section, the LTD is equipped in the open channel so that the relationship between the channel domain and the LTD domain is ∈ . (1) The LTD consists of fluid and tree domains, f and tree , as in which a Cartesian coordinate system Here, the subindex f represents the first letter of "fluid." Accordingly, the boundaries of the rectangular channel perpendicular to the ±x 1 -, ±x 2 -, and ±x 3 -directions are denoted as ±1 , ±2 , and ±3 , respectively. In this study, both air and water are assumed to be incompressible Newtonian fluids. The overall flow direction is also assumed to be in the x 1 -direction. For simplicity, trees are assumed to be rigid, that is, not deformable.
Although the deformation of some tree species might affect a lot in flow motion, 21 such rigid body assumption must be justified by referring previous laboratory experiments (e.g., Zhang et al. 33 ). These assumptions allow us to describe the fluid motion inside the channel by the following Navier-Stokes momentum and continuity equations: where is the fluid's mass density, u is the flow velocity vector, and f is the body force vector generally driven by the gravitational effect. Here, the stress tensor is expressed in terms of the pressure p and strain rate tensor d as where is the viscosity and I is the second-order identity tensor. The strain rate tensor d has been defined as the sum of the symmetric component of the velocity gradient as which is deviatoric due to incompressibility. Note that Equations (3)-(5) are applied not only to the flow in the LTD but also to the flow through the entire open channel domain . As mentioned in the previous section, the spatial arrangement of the LTD can be assumed to be periodic in the x 2 -direction such that the velocity field satisfies the periodic boundary condition in the transverse direction as In addition, the nonslip condition u = 0 is imposed on the bottom surface −3 and on the surfaces of the trees, denoted by tree .
On the upwind surface of the rectangular channel, the following inflow conditionû is applied on the surfaces perpendicular to the overall flow direction, that is, the x 1 -direction: where• indicates prescribed values. A Sommerfeld nonreflecting boundary condition 34 is applied to the outflow surface +1 by solving the following radiation equation where l indicates the time step of calculation. See, for example, Yoshida and Watanabe 35 for a detailed explanation of Sommerfeld nonreflecting boundary conditions. Fluid flow in the LTD is composed of water and air with an interface such that a two-phase flow is realized. Therefore, the fluid domain f is composed of the water flow domain w and air flow domain a as The following relationship is realized at the interface between the two phases.
where n is a unit normal vector orthogonal to the surface of each phase boundary and where w and a are boundaries of w and a , respectively, and are identified on the interface between the two phases. Here, the depth of water flow can be defined as where e +3 is the x 3 -direction unit normal vector. Thus, the inflow depthĥ, which is the vertical position of the interface, is prescribed to determine the flow rate of water in the inflow condition.
As explained previously, since the 3D NSs are carried out in the entire domain of the open channel , two inflow parameters, û andĥ, are applied to the boundary surface of the channel . However, only flow inside the LTD domain is utilized in the numerical flow test evaluation process. Thus, the inflow conditions for numerical flow tests in the LTD are not necessarily the same as those in the open channel.

Evaluation of the drag parameter at the global level
In this study, we propose a method for evaluating the global pressure loss of a coastal forest relevant to the global drag effect caused by trees in the LTD by conducting numerical flow tests (i.e., 3D NSs). In the hydraulic dynamic handbook written by Idelchik, 36 the pressure loss ΔP of the whole system due to certain resistances is generally defined as in a 1D setting, where  is the stream velocity and is the overall pressure loss coefficient that reflects the geometrical characteristics of the system.
In the 2D vegetated flow situation, with the assumption that the flow is unidirectional and the drag effect is isotropic from a macroscopic perspective, the global flow direction can be coincident with the x 1 streamwise direction without loss of generality. Thus,  in Equation (13) is identified with the x 1 -component of the global velocity vector, which is obtained by the volume-averaging operation for local variables over the LTD. In addition, the pressure loss in the x 1 -direction is highly influential on the drag effect compared with the loss in the other two directions. By referencing Figure 1C,D, Equation (13) can be used to determine the global pressure gradient P x 1 as where L = x out − x in is the distance between these two positions. We have also defined the drag parameter as C = ∕L. Additionally, the global pressure gradients in the x 2 -and x 3 -directions can be neglected in this study. Equation (14) would have validity considering that the mean pressure gradient is comparable to the sum of the bottom friction effects and the drag of trees under steady uniform flow conditions. 32 Thus, the goal in the third step is to compute the drag parameter C = ∕L by evaluating the global variables,  and P x 1 , from the numerical test results through the use of appropriate averaging operations, which will be explained as follows. It can be found that the parameter C is an alternative to the drag coefficients in the general drag force formulation. 37 Although the parameter C seems to work similarly to the classical drag coefficients, the parameter C can be interpreted as containing the information about the trees' density, their geometrical features and the dependence of the submerged depth as a result of numerical flow tests. We calculate the "global velocity"  by using the local velocity u distributed in the LTD as where we have defined the spatial and temporal averaging operations as Based on the theory of multiscale modeling, the global quantities calculated with these space-time-averaging operations in Equations (15) and (16) show the macroscopic flow characteristics of a homogeneous domain equivalent to the LTD with local heterogeneities. Here, | w | is the volume of water flow in the LTD, t 0 − T is the onset time for time integration, and T is a certain time interval. Following the idea proposed in previous research. 38 we specify T as Additionally, the following approximation can be made for the definition of | | considering the dominance of streamwise directional flow: Similarly, the global pressure gradient P x 1 is calculated as where we have defined the surface averaging operation as Here, +1 and −1 are the outflow and inflow surfaces of the LTD, respectively. By using the global velocity and pressure gradient evaluated with Equations (15) and (19), we can compute the drag parameter in Equation (14) as Furthermore, the average depth in the LTD, denoted by , can be calculated in a similar fashion as and is referred to as the "global depth" in this study.

Response surface as a surrogate
Since numerical flow tests are expensive and time consuming, only a small number of drag parameters C in Equation (21) can be obtained in response to various inflow conditions. However, each value of C calculated by a numerical flow test can be regarded as a function value defined by the velocity and depth that are representative of the flow conditions inside the coastal forest; this function is supposed to constitute a "response surface." Once this response surface is determined in some way, it can act as a surrogate for the costly numerical testing system that provides the drag parameter C as a function of these inflow parameters. This type of surrogate modeling has been utilized especially for engineering design. 39 In this study, an FEM-type approximation scheme is adopted to represent the function of the parameter C. Specifically, the space defined by the independent variables, which are arguments of the function, is discretized into a finite element (FE)-type mesh, which is composed of standard 2D elements. Then, each node corresponds to the set of arguments for a particular numerical flow test, and each element is a small area to interpolate the function with its nodal values. Thus, the function can be approximated by standard C 0 -continuous shape functions commonly used in FE approximations. The introduction of this kind of surrogate modeling must be one of the new contributions in this study. Let C * (s * ) be the response surface of the drag parameters under arbitrary inflow conditions s * = [ * ,  * ] and be approximated as where f (s i ) is the value of C calculated following the procedure in the previous subsection at the i-th control point representing a flow condition s i = [ , ], e is the subdomain of the whole parameter domain (= ⋃ e ), and N c is the number of control points in e . In this study, we employ triangular subdomains determined by three control points (N c = 3) to approximate the response surface so that linear interpolation functions can be utilized.

Inflow conditions
To demonstrate the proposed method, a representative numerical example is presented in this section. Although the whole process in this section is applicable for real trees, we conduct a series of numerical flow tests with a relatively small-scale LTD. A portion of the actual miniature forest is illustrated in Figure 2A, in which the staggered arrangement of miniature trees is located in the central part of the open channel. The domain surrounded by the red-colored solid line shown in Figure 2A is defined as the LTD for the numerical flow tests; see Figure 2B for a bird's-eye view with positions indicating the forest region and the boundaries. The LTD is established relatively remotely from both the inflow and the outflow surfaces −1 , +1 of the open channel. This is to avoid the undesirable influence of the inflow/outflow conditions on the flow inside the LTD. Figure 2C also shows the forest region, and Figure 2D illustrates the tree's geometric information, which is divided into three portions in view of the submergence depths. With the help of the stabilized FEM, 40,41 the 3D Navier-Stokes equation (3) is solved for numerical flow tests. The phase-field method 42,43 is also employed to capture the free surface of the flow. Precise explanations of these numerical schemes are presented in Appendix A. FE mesh discretization by the tetrahedron 4-node elements is illustrated in Figure 3 with the size of representative meshes summarized in Table 1. Here, the sizes of those meshes relative to the structure scale have been validated through both the previous studies. 44 Eight combinations of velocities and submerged levels are selected for the inflow conditions in the numerical flow tests within the open channel; their specific values are presented in Table 2. Two cases with higher (û 1 = 0.3 m/s) and lower (û 1 = 0.5 m/s) velocities are considered. When considering flow at a real scale, the Froude similarity law can be applied. Specifically, the inflow on trees at a real scale could be estimated by the application of the relationship between the Froude number and the inflow speed û 1 as In this study, all of the inflow conditions were set to realize "sub-critical flow" (Fr < 1). This setting is justified due to the fact that many of previous studies employed sub-critical flow situations to examine the tsunami flow through the vegetation. 45,46 x y x-y plane x-z plane However, as mentioned before, these levels are the input conditions for the open channel but not for the LTD that defines the region occupied by the coastal forest only. In this respect, the actual levels of the flow conditions inside the LTD will be defined in the next subsection.

Results of the numerical flow tests
A resultant flow in Case 1 is visualized in Figure 4. The array of woods is colored white, and the velocity distribution of the surface flow at each time step is provided. We can understand that the boundary conditions imposed on both the inflow and the outflow surfaces have affected the increased velocity at both edges, as shown in Figure 4A. We also understand that the flow, which is driven by the boundary condition at −1 , inundates the LTD zone from Figure 4B-D. After the flow passes through the trees at t = 2.0 s, the wake behind each tree can be observed, as shown in Figure 4E. The higher velocity observed behind the LTD at t = 2.0 s must be the result of wave encounters driven by both the inflow and the outflow surfaces. The flow speed increases locally in the space between downstream trees, as demonstrated in Figure 4E,F. We can also recognize that the flow speed is attenuated behind the branches. These observations indicate that the complex shape of branches can constitute a severe factor influencing the evaluation of flow attenuation. In contrast, the change during t = 10.0-20.0 s is less significant, as shown in Figure 4G-I. Figure 5 shows the 3D NS results in the LTD at t = 20.0 s obtained through a series of numerical flow tests with the inflow conditions summarized in Table 2. Here, the total time for each numerical flow test is set to 20.0 s. These  The initial flow depthĥ is established considering the submerged level corresponding to the vertical shape of trees as shown in Figure 2C. figures demonstrate that the 3D NSs with the stabilized FEM and phase-field method successfully capture the visible flow characteristics according to the inflow conditions. Indeed, these figures show that the nonuniformities of the velocity distributions in the cases with high inflow velocities (Cases 2, 4, 6, and 8) are more prominent than those with low inflow velocities (Cases 1, 3, 5, and 7). In particular, the results of Cases 2-8 are more turbulent than those of Case 1. This indicates that the cases with intermediate, higher, and maximum depths exhibit the influence of the canopy's complex geometry, which cannot be represented by 2D simulations or model experiments with arrays of simple cylinders.
Note that the flow depths and the velocities inside the forest region are different from those provided as inflow conditions in Table 2. For example, both Cases 1 and 2 were categorized into stem-submerged levels in Table 2, but the depth inside the forest obtained for Case 2 attains the canopy zone defined in Figure 2D as a result of the 3D NS; see Figure 5 for Case 2. Such a surface-increasing tendency is thought to be induced by the resistance effects of trees. As a result of the high flow velocity, the trees' resistance force is enhanced, thereby decreasing the permeability of the LTD. A similar tendency was observed in previous 3D simulations with a simple array of cylinders. 22 Let us calculate the global velocity  i and depth  as we defined in Equations (15) and (18). Here, temporal data of the space-averaged velocity |⟨u⟩| and [h] −3 are provided in Figures 6 and 7, respectively. It can be understood that both of those values dramatically vary for the first 10 s. On the other hand, their variations are relatively calm in the last 10 s (t = 10.0-20.0 s). Although both the velocity  i and depth  do not reach steady-state in a rigorous sense, the flow inside the LTD can be regarded as stable in the sense that the inertial effects is not dominant from the macroscopic perspective. Therefore, we conclude that we should carry out our time-averaging operation (16) with the data recorded during the last 10 s. We set t 0 in Equation (16)    With the above setting and by using the numerical test results as the input data for Equation (15), the global parameters are calculated, and the simulation cases modified in terms of the LTD are provided in Table 3. The global parameters  1 and  will be used as independent variables of the function representing a response surface of the drag parameter C, as introduced in Section 2.4 and approximated in Section 2.5.
With reference to the distinct flow conditions in Table 3, a further discussion might be possible from the numerical test results. As shown in Figure 4B of the cases with fully submerged depth, the flow structure in the lower zone (0 ≤ z ≤ 0.22) is quite different from that in the zone above it (0.22 < z). In fact, the flow velocity in the upper zone is higher than that inside the forest region. This tendency is in agreement with the actual observations in the hydraulic experiment. 47 The lower velocity inside the forest region must result from the local drag effect caused by tree canopies and trunks. It is thus expected that the flow field above the trees is minimally affected by their presence.
A similar tendency is applied to the comparison between the calculated flow structures in the canopy and trunk zones. Specifically, the flow velocity in the canopy zone is slower than that in the trunk zone. To confirm this inference, let us observe the local flow structure shown in detail in Figure 8, in which the streamlines around the single tree are located at (x 1 , x 2 ) = (1.04, 0.0) [m] in Cases 1 and 6. Note that other trees surrounding the visualized tree are arrayed as shown in Figure 8C, but are opaqued in Figure 8A,B for the focus on the single tree. The streamlines in the case   with a trunk-submerged depth ( Figure 8A) are not significantly disturbed by the trunk, while those in the case with a canopy-submerged depth ( Figure 8B) are disturbed considerably more by the branches and flow velocity in the trunk zone than are those in the canopy zone. Therefore, it appears reasonable to conclude that conducting numerical flow tests by 3D NSs in consideration of trees' geometrical features enables us to adequately capture complex flow structures that are relevant in evaluating the global drag effect of a coastal forest.

Response surface of the global drag parameter
With the simulation results of the numerical flow tests in the previous subsection, the global drag parameter C can be determined by using Equation (21), which is a function of the global velocity  and pressure gradient P x 1 . The global velocity was determined by using Equation (15) in the previous subsection; in this subsection, the global pressure gradients between the inflow and outflow surfaces of the LTD are calculated by the formula in Equation (19). For this purpose, we first provide in Figure 9 the time variations of the calculated global pressure losses ΔP(t) divided by the length L and adopt the same onset time t 0 and the same time interval T as before in Equation (16) to obtain the global pressure losses, which are the space-time-averaged values of the corresponding local pressure losses Δp(x, t). As shown in each graph in Figure 9, the overall pressure loss in the case of a higher inflow velocity is larger than that of a lower velocity. This dependence of the pressure loss on the flow velocity can also appear in the classical Morison formula for calculating the drag force caused by vertically standing structures in a flow field. In addition, all the graphs in Figure 9 exhibit a common tendency; that is, the variation in P x gradually becomes gentle and then reaches stable states after reaching the peak points, although there are some fluctuations.
In this context, the settings of T and t 0 for time averaging in Equation (17) might be validated by checking the global steadiness of the pressure gradient. For this purpose, we divide the global pressure losses ΔP(t) into time-averaged and time-fluctuating components, ΔP and ΔP ′ (t), respectively, as and then expect that

F I G U R E 9
Temporal variation in the local pressure gradients ΔP∕L(= P∕ x 1 ) between the inflow and outflow surfaces Here, ΔP ′ (t) can be calculated by subtracting ΔP from ΔP(t). The arithmetic mean of the eight calculations of ΔP ′ ranges from ±10 −16 to ±10 −14 Pa, which is sufficiently small compared to that of ΔP (ranging from −10 1 to −10 2 Pa). This negligible order of ΔP ′ guarantees that the flow field from t 0 to t 0 − T satisfies the global stable condition postulated in Equation (26). With the acquired guarantees of the space-time-averaging formula (Equation 16), we evaluated ΔP and calculated the values of the drag parameter C by using Equation (21), as presented in Table 3.
These values can produce the global flow characteristics caused by the presence of trees in the coastal forest. Nonetheless, the drag parameter C evaluated as a function of  1 and the pressure gradient P x 1 is difficult to use for global flow simulations realized by the 2D SW equations, whose independent variables are  1 and . In other words,  is more suitable than P x 1 as an argument of a function of the drag parameter C.
Based on the above discussions, the response surface of the drag parameter C can be suitably constructed as a function of  1 and . Figure 10 shows the surface thus constructed, which is approximated by linear interpolation functions as explained in Section 2.5. Note that seven nodal points are included in addition to the eight simulation results on the assumption that C = 0 1/m at  1 = 0 m/s, or/and  = 0 m. As shown in this figure, the drag parameter C, which represents a degree of resistance, tends to increase with increasing global depth and velocity. Additionally, the dependency of the drag parameter C on the global velocity appears to be monotonic, but that on the global depth is fairly complicated, as expected. From the contour plot shown in Figure 10A, specific evaluations of the dependency of C on the global depth and the global In conclusion, we have successfully determined the drag parameter C for the 2D SW equation by virtue of a 3D NS, which appropriately reflects the effect of the geometric information of (miniature) trees in a virtual coastal forest. These findings also suggest that the drag effect of trees can be expressed by only two variables (global velocity and global depth) without any other time-consuming parameterization. This argument is validated by 2D NSs realized by the SW equation enhanced by the constructed response surface or, equivalently, a "surrogate model" of the drag parameter C as addressed in the next section.

VERIFICATION ANALYSIS
A series of numerical analyses of 2D SW flows are carried out to verify the validity of the response surface of the global drag parameter C, which was obtained as a function of the global velocity and depth in the previous section.

Governing equations
An inundation flow involving vegetation is regarded as a global flow in this study and is assumed to be governed by 2D SW equations, which are known to be the most common and practical in the areas of coastal and tsunami engineering. The global domain of this 2D inundation flow is denoted bỹ∈  2 , which should be distinguished from that of the LTD, ∈  3 . A Cartesian coordinate systemx = (x 1 ,x 2 ) is introduced to describe the set of SW equations, which can be derived by depth-averaged operations on the 3D Navier-Stokes equations (Equation 3) along with some simplifications 48 as where H is the depth, g is the gravitational acceleration, and U i are the depth-averaged mean velocity components, each of whose magnitude is assumed to be identified with the global velocity. Additionally,  B i is the basal shear acting on the bottom surface and is assumed to have the following function form that is commonly employed for tsunami inundation analyses: 11 Here, n is the roughness coefficient depending on geographical features. Following the SW tsunami simulation with a drag force term, 2 S i in Equation (27), representing the trees' drag effect against the global flow, is assumed to depend on the global velocity as which can be equated with the global pressure gradients P∕ x i in view of Equation (14). In fact, C can be identified with the drag parameter introduced in Section 2 to characterize the drag force caused by trees in a coastal forest, and its function form was illustrated in the response surface at Figure 10. Note that Equation (29) is similar-but not identical-to the general drag formulation (e.g., Morison et al. 37 and Tanaka 49 ). A relatively simplified form is realized based on the drag parameter C, which comprehensively includes the information about the number of trees, the shape of trees, or so on. Note that, when there are no trees in the domain, drag parameter C in Equation (21) becomes zero as P x 1 must be approximately zero. As a result, the drag force term S i in Equation (29) disappears, so that Equation (27) becomes identical to the general 2D SW equation without the effect of trees. This was also confirmed from the result of an additional NS, although it is not provided here.

Reproduction of SW flow thorough a miniature tree model with a response surface
To verify the validity of the response surface (surrogate model) of the global drag parameter C in S i of Equation (27) that is implemented in the SW equations (Equation 27), 2D flow simulations are conducted with the stabilized FE scheme; see Appendix B for a detailed explanation of this scheme. The inflow and outflow conditions are established to be consistent with those in the 3D NSs performed in the previous sections so that the results can be compared with those of the 3D NSs. Figure 11 shows the 2D domain prepared for the SW simulations. The inflow and outflow surfaces, which, respectively, correspond to −1 and 1 in the LTD for the numerical flow tests in Section 3, are located at x 1 = −0.5 m and x 2 = 1.5 m, respectively. Moreover, the vegetation zone of the coastal forest is placed from x 1 = −0.0433 m to x 1 = 1.0825 m. Within this vegetation zone, the drag parameter C varies according to the global velocity and depth,  = [U]̃t ree and  = [H]̃t ree , respectively, both of which are averaged from the solutions U i and H of the SW equations (Equation 27), as the function form of C has been approximated by the response surface constructed in Section 3.3. Figure 12 shows the profiles of the obtained global free surfaces (or depths) and velocities around the vegetation zone with different inflow depths. Furthermore, the 3D NS results (Cases 1, 6, and 7) are presented for comparison purposes. Here, the profilesŪ(x),Ū(x), h(x), and H(x) have all been obtained by utilizing the same temporal averaging procedure as Equation (16), as explained in Section 2.4 (t 0 = 15.0s, T = L∕û 1 ).
Overall, the results of the SW simulations with the response surface (surrogate model) differ only slightly from the 3D NS results. In particular, as seen from Figure 12A,C, the global surface profiles of Cases 1 and 7 are in good agreement with the 3D NS results. Specifically, the depth reduction effect, which represents a crucial role of trees, is successfully captured by the 2D SW simulations owing to the surrogate model that we have constructed from the results of numerical flow tests. However, as shown in Figure 12B, the flow surface of Case 6 is lower than that of the 3D NS at the downstream edge.
In addition to the flow depth, the flow velocities are also investigated. The results obtained by the two schemes are compared on the right-hand side of Figure 12 and in Table 4, which summarizes the mean flow velocities inside the forest. These profiles illustrate that the global velocities obtained by the 2D SW flow simulations tend to increase as the depth decreases. This is because the flow fluxes on the inflow and outflow surfaces are the same in these cases. This tendency is especially noticeable in Cases 6 and 7 (see Figure 12B,C, respectively). Furthermore, similar to the tendency of the depth profiles, the flow velocity is overestimated at the downstream edge. Nonetheless, it is reasonable to conclude that each of the flow velocities obtained by the 2D SW simulations is in close agreement with the averaged value of the local flow velocity from the 3D NS as summarized in Table 4.

DISCUSSION
We believe a developed framework could evaluate the global drag effect of vegetation that has heretofore been considered mostly by empirical knowledge or simplified modeling in conventional or single-scale approaches. Our present approach resolves the problem remaining in those conventional modeling, such as the discrepancy between the depth arise and resistance reduction. Furthermore, a surrogate model based on 3D NSs could directly deal with the complex shape of natural trees, although most of the previous approaches require the two-dimensional conversion that lacking detailed morphology information.
As we previously demonstrated in Section 4, our proposed approach has an advantage in reasonable computational costs in the estimation of the streamwise flow profile. We have confirmed that the 2D SW simulations equipped with our surrogate model, in which the idea of multiscale modeling is partially implemented, are reliable to some extent in comparison with the 3D NSs. The 2D simulation with the proposed approach was completed in only several minutes with a single-core processor, although a trial 3D flow calculation required approximately 10 days to carry out in with whole coastal forests domain with parallel computations (1086 cores, CPU: 1.4 GHz/3.0464 TFlops).
Those advantages might be sustained for larger, in other words, real scale flow situations. The number of degrees of freedom (so-called DOF), which deeply influences the computational costs, is not increased if the scale becomes larger since the sufficient FE mesh size depends mainly on the ratio to the representative length of the structure. We also mentioned that the possibility of the miniature scale flow simulations as the alternative of larger-scale flow based on the Froude similarity law in Appendix C.
It is also to be noted that the present framework has some advantages in tree modeling compared to flume experiments. Indeed, they are expensive. In addition, the realistic "physical" tree modeling in laboratory, which may require some handcrafting, is still time-consuming, even though the recent development of 3D-scanning and printing technologies simplifies the experimental setup. If we rely on a more sophisticated algorithm (e.g., L-system 50 ), our tree modeling procedure might be completed only in the workspace of hardware. Nevertheless, a real-scale experiment is unrealistic. On the other hand, the LTD in our numerical experiments is scalable to a realistically larger domain as the periodic boundary condition are imposed while the channel in a flume experiment is closed.
In the meantime, the proposed approach has some room for improvement. In particular, examining the 2D velocity profiles in detail in Section 4, we find that the flow velocities in the cases with relatively high depths (Cases 6 and 7  underestimated in comparison with the corresponding 3D NS results. This remains to be resolved in future studies by taking into accounts the influence of turbulence in more detail. Comparison with some laboratory experiments is also needed to demonstrate the capability of the proposed method. There is also a limitation such that the current framework could only deal with the stable, quasi-steady and macroscopically unidirectional flow situation. We formulate SW equation (27) without inertial force effects like many previous studies. 45 However, the inertial force is dominant in dense forests under unsteady-state flow conditions especially at the tsunami runup. 51 Since the flow has been assumed to be macroscopically unidirectional in this study, we have to investigate whether or not this assumption would cause a discrepancy by applying more complex inflow conditions in the macroscale SW simulations. This is also left to future work.

CONCLUSION
The main purpose of this study is to develop a method for introducing the resistance effect of trees in a coastal forest into global SW simulations by using the results of 3D NSs on the local level. Along the lines of multiscale modeling, the 3D NSs are referred to as "numerical flow tests" and are carried out in a LTD involving a coastal forest with various flow conditions to evaluate the 3D local flow characteristics, which are consolidated into the trees' drag performance in the global flow with the proposed volume-time-averaging procedure. Specifically, the pressure loss measured over the LTD volume and a certain time interval is regarded as the global resistance effect of the group of trees and then, with the help of a classic hydraulic formula, converted into the drag parameter C, which represents a resistance force in the global 2D SW equations. With the help of surrogate modeling, the drag parameter C is represented by a function or, equivalently, a response surface consisting of two flow variables: the global velocity  and global depth . Since the number of numerical flow tests is limited, the response surface has to be approximated from a discrete data set consisting of C,  , and  to be equivalent to the actual drag effect without any other time-consuming parameterization.
The procedure and capability of the proposed method were demonstrated by a numerical example. The obtained response surface, which was an interpolated function of the global drag parameter C, could reflect the complex flow characteristics of 3D flows involving trees in the LTD and was successfully utilized in the 2D SW simulations. To be more specific, the profiles of the global flow surfaces (depths) and the global velocities obtained by the 2D SW simulations were in close accordance with the 3D NS results, although slight underestimations were observed, especially when the depth and velocity increased.
Our present framework resolved the discrepancy remaining in much previous vegetation modeling. Our 3D NSs based framework could directly evaluate the complex shape with remarkable computational efficiency in 2D SW. Although currently there are some tasks and limitations that we have to resolve for making our approach to become more robust and practical, there is also the potential that it would contribute reasonable and reliable flow simulations involving coastal forests in the future.

ACKNOWLEDGMENTS
This work is supported by Grant-in-Aid for JSPS Research Fellow (JP17J05235) and JSPS KAKENHI (JP25246043). Also, the computations in this work have been done using the Fujitsu PRIMERGY CX1640M1 (Oakforest-PACS) at the University of Tokyo and the Cray CS400 2820XT (Laurel 2) at Kyoto University.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author, R. Nomura, upon reasonable request.
The fourth term is the SUPG/PSPG term added for stabilization. Here, supg and pspg are both stabilization parameters for the momentum equation obtained by the following formulation: The weak form of Equation (A9) is given in the following form with the test function * : where e denotes the domain of each element after FE discretization. Based on the SUPG scheme, the stabilization parameter is defined as For Equations (A2) and (A13), we employ the Crank-Nicolson implicit scheme because of its numerical stability. P1/P1 type elements are used in the spatial discretization.

B.1 Stabilized FEM for 2D SW Flow Simulations
The formulation for the 2D SW equation based on the stabilized SUPG scheme is derived in the same way as for the compressible flows provided by Tezduyar et al.; 52 see Takase et al. 53 for the previous application of this formulation. With the state vector U = [H, U 1 H, U 2 H] T , the governing equations (27) are rewritten in the following form: where A i (i = 1, 2) are flux Jacobian matrices given as follows: Here, c is the velocity of the gravity wave √ gH. The diffusion coefficient matrices K ij , (i, j = 1, 2) are defined as The eddy viscosity is calculated by the following form: Here, is the von Karman constant, for which the generally accepted value is 0.41, and u * is the friction velocity, generally defined with the Manning roughness coefficient n: R and G are the stress tensors caused by pressure and bottom friction, respectively, and are given as The definition of  B i follows Equation (28). The porous drag term S is given as follows: After substituting the test function U * into Equation (B1), the resulting weak form can be obtained: Here, SW supg and SW const are stabilization parameters defined as follows: In addition, ||u h ||, k , and zare calculated at each time step as follows: (B14)

APPENDIX C. APPLICABILITY TO REAL SCALE FLOW
Although this paper demonstrated only by the 0.22 m height miniature tree model shown in Figure 2D, our numerical test would be valid for larger identical scale models by using Froude similarity law. Here, relationship between the flow through miniature model and real scale model can be described by the following formulation: where û 1,mini and û 1,real means the inflow speed controlled in numerical flow test and equivalent flow speed in times larger real scale domain, respectively. Here, inflow depthĥ mini andĥ real could be simply related tô According to Equations (C1) and (C2), a variety of real flow parameter could be mutually translated between real scale and miniature scale.
For instance, we setup the isolated tree in both the channel described in Figure 2B and 10 times larger scale channel domain ( = 10) as shown in Figure C1. While the inflow conditions of case at Table 2 can be applied for miniature channel with isolated models. With the relationship stated in Equations (C1) and (C2), the equivalent inflow and other numerical input conditions for 10-times larger scale model can be identified in Table C1.
The results obtained in both two scale model are provided in Figures C2 and C3. The vertical flow profiles behind the isolated tree model observed in two scales are provided in Figure C2. As can be seen from the vertical profile behind the isolated tree model Figure C2 Figure C2A into the equivalent real scale flow by Froude similarity as shown in the solid line in Figure C2C. We could clearly understand that the tendency is quite similar between the two simulation results.

F I G U R E C1
Real scale channel domain (10 times larger = 10 than miniature channel in Figure 2D) with an isolated tree model In addition to that, we calculated the fluid forces exerted on rigid trees as the quantitative evaluation indices from the following equation: The evaluated values from Equation (C3) are illustrated in Figure C3A,B. Based on the Froude similarity, the force exerted in times larger scale is equivalent to 3 values measured in its miniature scale. This 3 relationship seems to be satisfied when we compare two graphs. Such insights can be further confirmed by plotting two data at the same graphs. The dashed line in Figure C3C shows the 3 times force values at 1∕2 larger time scale provided in Figure C3A. The dashed line follows the solid lines, which indicates the value measured in a real scale model ( times large model). Although slight differences appear as the time elapsed, both model shows similar trend especially in the maximum peak values around t = 10 s. Such insights are further convinced by referring to the time averaged values of force measured in both two scales, summarized in Table C2. The force values exerted on real scale model simulations are almost equivalent to that estimated from 3 relationship deployed by Froude similarity law.