A comparison of a two-dimensional depth-averaged flow model and a three-dimensional RANS model for predicting tsunami inundation and fluid forces

The numerical modeling of tsunami inundation that incorporates the built environment of coastal communities is challenging for both 2-D and 3-D depth-integrated models, not only in modeling the flow but also in predicting forces on coastal structures. For depth-integrated 2-D models, inundation and flooding in this region can be very complex with variation in the vertical direction caused by wave breaking on shore and interactions with the built environment, and the model may not be able to produce enough detail. For 3-D models, a very fine mesh is required to properly capture the physics, dramatically increasing the computational cost and rendering impractical the modeling of some problems. In this paper, comparisons are made between GeoClaw, a depth-integrated 2-D model based on the nonlinear shallow-water equations (NSWEs), and OpenFOAM, a 3-D model based on Reynolds-averaged Navier–Stokes (RANS) equation for tsunami inundation modeling. The two models were first validated against existing experimental data of a bore impinging onto a single square column. Then they were used to simulate tsunami inundation of a physical model of Seaside, Oregon. The resulting flow parameters from the models are compared and discussed, and these results are used to extrapolate tsunami-induced force predictions. It was found that the 2-D model did not accurately capture the important details of the flow near initial impact due to the transiency and large vertical variation of the flow. Tuning the drag coefficient of the 2-D model worked well to predict tsunami forces on structures in simple cases, but this approach was not always reliable in complicated cases. The 3-D model was able to capture transient characteristic of the flow, but at a much higher computational cost; it was found this cost can be alleviated by subdividing the region into reasonably sized subdomains without loss of accuracy in critical regions.

1 Introduction

1.1 Scope and motivation of the study

For many years, researchers have been working on different numerical models that can predict tsunami behavior. Tsunami prediction generally requires modeling at a wide range of spatial scales, including (from large to small scale) offshore wave propagation, beach run-up, inland inundation, and impact on individual structures.

Due to the large differences in scale for the different processes, most tsunami models solve two-dimensional depth-integrated equations, e.g., the nonlinear shallow-water equations (NSWEs) or some form of Boussinesq wave equations, to predict tsunami behavior, using computational grids that vary several orders of magnitude in spatial resolution, from several kilometers far from the shoreline to 10 m inland. The NSWEs are often used in the nearshore and inundation zone, since they can handle nonlinearities that arise in very shallow water and can be adapted to deal robustly with wetting and drying. However, it is not clear whether these equations are adequate to properly model fully three-dimensional turbulent flow, particularly at the scale necessary to determine tsunami impact and corresponding tsunami-induced forces on individual structures.

It would be preferable to solve the three-dimensional Navier–Stokes equations with a proper turbulence closure. However, this is still extremely expensive computationally relative to two-dimensional models, and only practical for detailed simulations over small spatial regions.

Although such modeling is challenging, the latest version of ASCE 7-16 (American Society of Civil Engineers) in the United States has a chapter for tsunami loads and effect for coastal structures. The provision requires site-specific inundation modeling and analysis be performed for all vertical evacuation structures. One of such examples is the design of the first vertical evacuation structure in the United States (Ash, 2015), the site-specific inundation analysis of which was conducted by González et al. (2013).

Given the challenges in tsunami inundation modeling, and the necessity of doing so for coastal structures design and tsunami hazard mitigation, the goal of this paper is to (1) study the characteristics of two different types of models when applied to tsunami inundation modeling (with an explicitly represented constructed environment); (2) show some guidance for modeling tsunamis or other flooding events in similar constructed environments; and, by showing the pros and cons of both models, (3) provide some insight to tsunami hazard researchers, coastal structure designers, and relevant government agencies.

1.2 Previous work

Before introducing the two numerical models used in the current study, a brief review of previous research involving different types of models is given below.

The two-dimensional depth-integrated equations are the most widely used tsunami models for their simplicity and computational efficiency. Popinet (2012) simulated the 2011 Tōhoku tsunami by solving the 2-D NSWE with dynamically adapted spatial resolution that varied from 250 m in flooded areas nearshore up to 250 km offshore. The model accurately predicted long-distance wave and coarse-scale flooding; the initial surface elevation was determined from a source model based on seismic inversion (as opposed to inversion of DART buoys and tidal gauge time series). This also showed that an accurate and consistent model of tsunami wave propagation can sometimes be constructed using only seismic wave inversion. Wei et al. (2013) used the Method of Splitting Tsunamis (MOST) model to simulate the same tsunami event. The MOST model solves the shallow-water equations in spherical coordinates with numerical dispersion. Their results demonstrated that it may be possible to forecast near-field tsunami inundation in real time. Hu et al. (2000) presented an NSWE model that can simulate storm waves propagating in the coastal surf zone and overtopping a sea wall. They found that waves overtopping a vertical wall may be approximately modeled by representing the wall as a steep slope and that the overtopping rate is sensitive to the bottom friction and the minimum friction depth. The two-dimensional NSWE model of wave run-up and overtopping by Hubbard and Dodd (2002) features an adaptive mesh refinement (AMR) algorithm. Their model can accurately reproduce 1-D and 2-D wave transformation, run-up, and overtopping in physical experiments. Their modeling of seawall overtopping by off-normal incident waves showed that there can be more flooding in such a situation than at normal incidence. Lynett (2007) simulated long-wave run-up obstructed by an obstacle and concluded that the obstacle can help reduce run-up and maximum overland velocity if the wave is highly nonlinear (with a ratio of wave height to shelf water depth ≥0.5). The sensitivity study also showed that in cases of breaking waves the Boussinesq model was more accurate than the nonlinear shallow-water equations in terms of wave run-up (maximum differences up to 10 %). For nonbreaking long waves, differences between the two were negligible. Shi et al. (2012) developed a high-order adaptive time-stepping total variation diminishing (TVD) solver for a fully nonlinear Boussinesq model and validated it against a series of laboratory experiments for wave shoaling and breaking and a suite of benchmark tests for wave run-up. The results showed that the model was able to accurately model wave shoaling, wave breaking, and wave-induced nearshore circulation. With a Boussinesq model, Lynett et al. (2010) simulated overtopping of levees of the Mississippi River–Gulf Outlet (MRGO) during Hurricane Katrina at four characteristic transects along the 20 km long stretch of the levees. The predicted overtopping rates agreed well with the observed data.

As computing power increases, it becomes possible to model the tsunami run-up process by solving three-dimensional Navier–Stokes equations with a proper turbulence closure. Choi et al. (2007) solved three-dimensional Reynolds-averaged Navier–Stokes (RANS) equations to simulate wave run-up on an conical island and compared different turbulence closure models including k-ϵ, RNG (renormalization group methods; Yakhot et al., 1992), and LES (large-eddy simulation). Their results showed that LES and RNG k-ϵ are similar and more accurate than k-ϵ is worse than those two. Williams and Fuhrman (2016) solved incompressible RANS equations with a transitional variant of the standard two-equation k-ω turbulence closure to study boundary layer flow induced by tsunami-scale waves. Their results indicated that the boundary layer generated by a tsunami is both current-like due to the long duration and wave-like due to its unsteadiness. The study also indicated that an existing expression for maximum bed shear stress under wind wave scale can be reasonably extrapolated to full tsunami scale. Mayer and Madsen (2000) investigated wave breaking in the surf zone by solving the RANS equations with a k-ω turbulence model. They found that the volume-of-fluid method could be used successfully to simulate wave breaking and that, although some instabilities occurred in applying the RANS equations, they can be eliminated by an ad hoc modification of the turbulence model. Other relevant studies include Biscarini (2010), Montagna et al. (2011), and Larsen et al. (2017), which model landslide generated tsunamis and tsunami-induced scours around coastal structures.

The prediction of tsunami impact on individual structures is also important because it provides guidance on designing coastal structures in tsunami inundation zones. The two-dimensional depth-integrated model may not work properly for these scenarios since the problems are more three-dimensional with large variation in the vertical direction and with transient and turbulent flow impacting the structure. In these cases, a three-dimensional model that solves the Navier–Stokes equation may give much better results. Researchers at the University of Washington (UW) modeled a series of “dam break” experiments by solving the 3-D RANS equations for bore-type impact of a wave on a series of 1∕20-scale model girder bridges to assess the 3-D effects on bridge skew (Motley et al., 2015; Wong, 2015).

The scale of modeling tsunami inundation inland with an explicitly represented constructed environment lies between that of modeling the large-scale tsunami wave propagation offshore and the small-scale tsunami impact on individual structures. This process is actually even more challenging to model since, for two-dimensional depth-integrated models, inclusion of the constructed environment increases the complexity of the topography, and the flow begins to have more variation in the vertical direction; while for the three-dimensional model that solves the Navier–Stokes equations, a fine mesh needs to be generated around each individual structure, which dramatically increases the number of cells in the computational domain.

Some researchers have tried to model this process with two-dimensional models. Ozer Sozdinler et al. (2015) used the numerical code NAMI DANCE to investigate hydrodynamic parameters in tsunami inundation zones with idealized structures – three rows of 20 blocks representing three-story concrete buildings. The code solved the NSWE using a finite-difference technique in a staggered leapfrog scheme. The effect of wave period, wave shape, protection structures, building layout, and Manning's coefficient are discussed. Some major conclusions included that the coastal protection structures like seawalls and breakwaters have very limited effect if the waves are able to overtop them and that it is preferable to use different Manning's coefficients for the sea, land, and buildings if more accurate values of hydrodynamic parameters are needed, albeit at the expense of more computational time. Similar conclusions on the Manning's coefficient were presented by Park et al. (2013). They simulated tsunami inundation in part of Seaside, Oregon, and compared flow parameters with their physical experiment. The comparison showed that the flow parameters were sensitive to the friction coefficient, especially for the momentum flux, which is proportional to tsunami loads on structures. For instance, decreasing the friction coefficient by a factor of 10 increased the predicted momentum flux by 208 %. Muhari et al. (2011) compared three different tsunami inundation models to evaluate tsunami impact on coastal communities: (1) the Constant Roughness Model (CRM) which uses a constant friction coefficient, does not include the constructed environment, and assumes that all buildings are not able to withstand the tsunami; (2) the Topographic Model (TM), which includes the constructed environment by incorporating building shape and height information into the topography; and (3) The Equivalent Roughness Model (ERM), which represents the building by using a different equivalent friction coefficient at the site of a building on the original topography (with only terrain information but not building height). Both the TM model and the ERM model gave more reliable prediction than the CRM model did, which confirmed the importance of taking the constructed environment into consideration.

However, few researchers have tried to use a three-dimensional model for inundation in a complex built environment. Shin et al. (2012) applied a 3-D LES model with two-phase flow to simulate inland tsunami inundation in a coastal city with hundreds of buildings and compared the prediction with experimental measurements. However, a fairly coarse mesh was used on land, and each building had only three to five mesh cells along its edge in the along-shore or cross-shore direction, so that the resulting agreement in flooding depth can only be considered qualitative. Qin et al. (2018) used 3-D RANS equations to predict tsunami inundation process and loads on individual buildings in part of Seaside, Oregon, and demonstrated that the whole part can be modeled using subsections with proper width without loss of accuracy in areas of interest.

In this paper, the results from a two-dimensional NSWE model and a 3-D Navier–Stokes model are presented for the test case of flow through a scale model of a portion of Seaside, Oregon. Two open-source models are used: the 2-D GeoClaw software from Clawpack (Clawpack Development Team, 2015), which is widely used for modeling tsunamis (both global propagation and local inundation), and the 3-D OpenFOAM software (The OpenFOAM Foundation, 2014). The two models are first compared and validated against an experiment in which a simple bore impinges on a single column and then compared for the Seaside model.

2 Simulation methodology

2.1 Two-dimensional model

The nonlinear shallow-water equations can be written as

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{h}_t+(uh{)}_x+(vh{)}_y=\mathrm{0},\text{(2)}& {\displaystyle }& {\displaystyle }(hu{)}_t+(huv{)}_y+{\left(h{u}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}\right)}_x=-gh{B}_x-Du,\text{(3)}& {\displaystyle }& {\displaystyle }(hv{)}_t+(huv{)}_x+{\left(h{v}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}\right)}_y=-gh{B}_y-Dv,\end{array}$$

where $u(x,y,t)$ and $v(x,y,t)$ are the depth-averaged velocities in the two horizontal directions, h is the water depth, g is gravitational acceleration, B(x,y) is elevation of the seabed, and $D=D(h,u,v)$ is the drag coefficient. The drag coefficient D could have many forms; in this study it is represented by

$$\begin{array}{}\text{(4)}& D={\displaystyle \frac{g{M}^{\mathrm{2}}\sqrt{(u^{\mathrm{2}}+v^{\mathrm{2}})}}{h^{\mathrm{1}/\mathrm{3}}}},\end{array}$$

where M is the Manning's coefficient and is set to 0.025 s m${}^{-\mathrm{1}/\mathrm{3}}$ for all two-dimensional simulations in this study. The subscripts in these equations represent first-order partial derivatives.

The GeoClaw model (LeVeque et al., 2011; Berger et al., 2011) features AMR and is released as a submodule of the Clawpack software (Clawpack Development Team, 2015), an open-source package for solving hyperbolic systems of partial differential equations (PDEs) of one, two, and three dimensions, through finite-volume implementation of high-resolution Godunov-type “wave-propagation algorithms”. Cell averages of the solution variables q are computed over the volume of each cell and updated with waves propagating into the cell from all surrounding cell edges. The wave at each edge is computed by solving a “Riemann problem” with initial piecewise constant data determined by cell averages on each side of the edge. This method is especially good at solving problems with discontinuous solutions like shock waves, which usually arise in the solution of nonlinear hyperbolic equations (e.g., bores in the case of NSWEs).

Specifically, GeoClaw uses a variant of the f-wave formulation of the wave-propagation algorithms that allow incorporation of the topography source terms on the right-hand side of Eqs. (2) and (3) into the Riemann problem directly. The augmented Riemann solver in GeoClaw combines the desirable qualities of the Roe solver (Roe, 1981), HLLE-type (Harten, Lax, van Leer and Einfeldt) solvers (Einfeldt, 1988; Einfeldt et al., 1991), and the f-wave approach (Bale et al., 2003). The Roe solver provides an exact solution for the single-shock Riemann problem. It is also depth-positive semidefinite like the HLLE solvers, has a natural entropy fix by providing more than two waves, and yields a better approximation for problems with large rarefactions. A large class of steady states is also preserved, even for non-stationary steady states with nonzero fluid velocity. In addition, it is able to handle the presence of dry states in the “Riemann problem”, in which one state is wet (h>0) while another is dry (h=0), or both states are dry. It also works robustly in situations where the topography changes abruptly from one cell to another by an arbitrarily large value. For more details on the augmented Riemann solver in GeoClaw, see George (2008).

A typical characteristic of tsunami inundation models, especially those that incorporate the built environment, is that the spatial scale of regions of interest may vary from kilometers to meters. For regions several kilometers offshore, grid cells can be as large as thousands of meters on a side, while for regions near the shoreline or in a built environment onshore, grid cells must be refined to several meters or less, since the size of a building may be only several meters and an adequate number of grid cells are required to achieve acceptable accuracy. In GeoClaw, a patch-based AMR technique can efficiently handle these situations (LeVeque et al., 2011; Berger and Leveque, 1998).

2.2 Three-dimensional model

For the three-dimensional model, version 2.3.1 of the open-source computational fluid dynamics (CFD) package OpenFOAM was used (The OpenFOAM Foundation, 2014). The package comes with different solvers for different types of flow. For tsunami inundation, in which there are two immiscible fluids (air and water) with a free interface, the interFoam solver can be chosen which uses the PISO algorithm to solve the RANS equations with a volume-of-fluid (VOF) approach to model the free surface. For details on these numerical methods, readers can refer to Hirt and Nichols (1981) and Versteeg and Malalasekera (2011). The VOF approach defines a scalar field α_water, which represents fractional volume of water in each cell. A cell full of water (ρ=1000 kg m⁻³, $\mathit{\nu }=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹) has α_water=1.0, while a cell full of air (ρ=1.22 kg m⁻³, $\mathit{\nu }=\mathrm{1.48}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹) has α_water=0.0. Here ρ is the mass density of the fluid and ν is the kinematic viscosity. A cell with α_water between 0 and 1 contains the interface. A special transport equation is solved to advance the α_water field. To close the RANS equations, Menter's k-ω shear stress transport (SST) model (Menter and Esch, 2001) was applied.

There are many other turbulence closure models, among which the k-ϵ model is also very popular. It is suitable for fully turbulent and non-separated flows and has the shortcoming of numerical stiffness in the viscous sublayer, which can result in stability issues (Menter, 1993). It was also applied to model the inundation process in this study but became unstable during the simulation. The k-ω SST is generally more stable and behaves better in modeling partially separated flows, which is the case in the current study (flow becomes separated after passing around the built environment).

With the assumption of an incompressible fluid, the RANS equations are listed below:

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial x_i}}=\mathrm{0},\text{(6)}& {\displaystyle }& {\displaystyle }\mathit{\rho }{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial t}}+\mathit{\rho }{\stackrel{\mathrm{‾}}{u}}_j{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial x_j}}=-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial x_i}}+\mathit{\mu }{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u}}_i}{\partial x_j\partial x_j}}-{\displaystyle \frac{\partial \mathit{\rho }\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}}{\partial x_j}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{u}}_i$ is the mean velocity in the i direction, ${u_i}^{\prime }$ is the fluctuating component of velocity in the i direction, and $\stackrel{\mathrm{‾}}{p}$ is the mean pressure. If u_i is the velocity component in the i direction, then $u_i={\stackrel{\mathrm{‾}}{u}}_i+{u_i}^{\prime }$. The Reynolds stress term in Eq. (6) is

$$\begin{array}{}\text{(7)}& -\mathit{\rho }\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}={\mathit{\nu }}_{\mathrm{t}}\mathit{\rho }\left[{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_j}{\partial x_i}}\right]-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}k\mathit{\rho }{\mathit{\delta }}_{ij},\end{array}$$

where k is the turbulence kinetic energy and ν_t is the turbulence eddy viscosity. The equations above need to be closed with some closure model. Here Menter's k-ω SST model (Menter and Esch, 2001) was applied:

$$\begin{array}{ll}\text{(8)}& {\displaystyle }& {\displaystyle \frac{\partial k}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathbf{U}k\right)=\stackrel{\mathrm{̃}}{G}-{\mathit{\beta }}^{\ast }k\mathit{\omega }+\mathrm{\nabla }\cdot \left[\left(\mathit{\nu }+{\mathit{\alpha }}_k{\mathit{\nu }}_{\mathrm{t}}\right)\mathrm{\nabla }k\right],{\displaystyle }& {\displaystyle \frac{\partial \mathit{\omega }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathbf{U}\mathit{\omega }\right)=\mathit{\gamma }{S}^{\mathrm{2}}-\mathit{\beta }{\mathit{\omega }}^{\mathrm{2}}+\mathrm{\nabla }\cdot \left[\left(\mathit{\nu }+{\mathit{\alpha }}_{\mathit{\omega }}{\mathit{\nu }}_{\mathrm{t}}\right)\mathrm{\nabla }\mathit{\omega }\right]\\ \text{(9)}& {\displaystyle }& {\displaystyle }+\left(\mathrm{1}-F_{\mathrm{1}}\right){\mathrm{CD}}_{k\mathit{\omega }},\end{array}$$

where ν is the kinematic viscosity of fluid and $\stackrel{\mathrm{̃}}{G}$ is defined as $\stackrel{\mathrm{̃}}{G}=min\left\{G,c_{\mathrm{1}}{\mathit{\beta }}^{\ast }k\mathit{\omega }\right\}$, where G is the production term and defined as

$$\begin{array}{}\text{(10)}& G={\mathit{\nu }}_{\mathrm{t}}{S}^{\mathrm{2}};\end{array}$$

S is the invariant measure of the strain rate, defined by

$$\begin{array}{}\text{(11)}& S=\sqrt{\mathrm{2}{S}_{ij}{S}_{ij}};\end{array}$$

and S_ij is the strain rate tensor, defined by $S_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\nabla }\mathbf{U}+{\mathbf{U}}^T\right)$. F₁ is a blending function defined by

$$\begin{array}{}\text{(12)}& F_{\mathrm{1}}=\tanh \left\{{\left\{min\left[max\left({\displaystyle \frac{\sqrt{k}}{{\mathit{\beta }}^{\ast }\mathit{\omega }y}},{\displaystyle \frac{\mathrm{500}\mathit{\nu }}{y^{\mathrm{2}}\mathit{\omega }}}\right),{\displaystyle \frac{\mathrm{4}{\mathit{\alpha }}_{\mathit{\omega }\mathrm{2}}k}{{\mathrm{CD}}_{k\mathit{\omega }}^{\ast }{y}^{\mathrm{2}}}}\right]\right\}}^{\mathrm{4}}\right\},\end{array}$$

where ${\mathrm{CD}}_{k\mathit{\omega }}^{\ast }$ is defined by

$$\begin{array}{}\text{(13)}& {\mathrm{CD}}_{k\mathit{\omega }}^{\ast }=max\left({\mathrm{CD}}_{k\mathit{\omega }},{\mathrm{10}}^{-\mathrm{10}}\right)\end{array}$$

and CD_kω is defined by

$$\begin{array}{}\text{(14)}& {\mathrm{CD}}_{k\mathit{\omega }}=\mathrm{2}{\mathit{\sigma }}_{\mathit{\omega }\mathrm{2}}\mathrm{\nabla }k\cdot {\displaystyle \frac{\mathrm{\nabla }\mathit{\omega }}{\mathit{\omega }}}.\end{array}$$

After solving Eqs. (8) and (9), ν_t can be calculated by

$$\begin{array}{}\text{(15)}& {\mathit{\nu }}_{\mathrm{t}}={\displaystyle \frac{a_{\mathrm{1}}k}{max\left(a_{\mathrm{1}}\mathit{\omega },S{F}_{\mathrm{2}}\right)}},\end{array}$$

where F₂ is a second blending function, defined as

$$\begin{array}{}\text{(16)}& F_{\mathrm{2}}=\tanh \left\{{\left[max\left({\displaystyle \frac{\mathrm{2}\sqrt{k}}{{\mathit{\beta }}^{\ast }\mathit{\omega }y}},{\displaystyle \frac{\mathrm{500}\mathit{\nu }}{y^{\mathrm{2}}\mathit{\omega }}}\right)\right]}^{\mathrm{2}}\right\}.\end{array}$$

All other constants are computed using a blend from the corresponding constants associated with the k-ϵ and k-ω models via blending functions like $\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{1}}{F}_{\mathrm{1}}+{\mathit{\varphi }}_{\mathrm{2}}\left(\mathrm{1}-F_{\mathrm{1}}\right)$. Values for these constants are α_k1=0.85013, α_k2=1.0, α_ω1=0.5, α_ω2=0.85616, β₁=0.075, β₂=0.0828, γ₁=0.5532, γ₂=0.4403, ${\mathit{\beta }}^{\ast }=\mathrm{0.09}$, a₁=0.31, and c₁=10.0 (Menter et al., 2003).

A force vector, F, on a structure is computed by summing forces from pressure, F_p, and from viscous stress, F_v.

$$\begin{array}{}\text{(17)}& \mathbf{F}={\mathbf{F}}_{\mathrm{p}}+{\mathbf{F}}_{\mathrm{v}}\end{array}$$

F_p and F_v are calculated respectively by

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }{\mathbf{F}}_{\mathrm{p}}=\sum _i\left(-p_i{A}_i{\left({\mathit{\alpha }}_{\mathrm{water}}\right)}_i{\mathit{n}}_i\right),\text{(19)}& {\displaystyle }& {\displaystyle }{\mathbf{F}}_{\mathrm{v}}=\sum _i\left\{\left({\mathit{\tau }}_i\cdot {\mathit{n}}_i\right)A_i{\left({\mathit{\alpha }}_{\mathrm{water}}\right)}_i\right\},\end{array}$$

where i is the index of cell faces on the building on which forces need to be evaluated; p_i is the total pressure on face i; A_i is area of face i; ${\left({\mathit{\alpha }}_{\mathrm{water}}\right)}_i$ is volume fraction of water in the adjacent cell of face i; n_i is the unit normal vector of face i pointing into the computational domain; and τ_i is the viscous stress tensor at face i, which can be expressed by ${\mathit{\tau }}_i=\left\{\mathit{\rho }\left(\mathit{\nu }+{\mathit{\nu }}_{\mathrm{t}}\right)\left[\mathrm{\nabla }\mathbf{U}+\mathrm{\nabla }{\mathbf{U}}^T\right]\right\}$ on face i.

3 Initial comparison of the 2-D and 3-D numerical models

An initial comparison of the two numerical models was conducted by modeling the interaction between a bore and a free-standing coastal structure, with experimental results from Árnason (2005). The experiment was performed at the Charles W. Harris Hydraulics Laboratory at UW, Seattle. In the experiment, a square column was placed in a 16.6 m long, 0.6 m wide, and 0.45 m deep wave tank, and aligned in parallel to the tank side walls (Fig. 1).

A thin gate separated water in the tank into two parts with different depths: 0.02 m deep on the square column side and 0.25 m deep on the other side. When the gate was lifted to the top of the tank in 0.2 s by a 6.4 cm diameter pneumatic piston, a bore formed and propagated toward the square column downstream. The square column with a 12×12 cm square-shaped cross section was placed 5.2 m downstream from the gate. To measure hydrodynamic forces, the column was supported from above and connected with a force sensor.

Both the three-dimensional and two-dimensional models were developed at model scale to simulate the physical experiment. The three-dimensional OpenFOAM model incorporated the column into the computational domain by simply cutting off a block of mesh of the same shape from the computational domain. The mesh was coarse far from the column (1 cm by 1 cm by 0.5 cm in the x, y, and z directions, where z is the vertical direction) and was refined gradually to 0.125 cm by 0.125 cm by 0.0625 cm in the x, y, and z directions near the column surface. These distances are evaluated to be 70, 70, and 35, respectively, in terms of dimensionless wall distance defined as $y^{+}=\frac{y{u}^{\ast }}{\mathit{\nu }}$, where y is the distance, u^∗ is the friction velocity, and ν is the kinematic fluid viscosity. The mesh was finer in the z directions to better capture the water surface. Forces on the column were obtained by integrating pressure and shear forces from fluid on the surface of the column.

Figure 1 Schematic of the experimental setup for the interaction between bore and square column. The top figure shows a plan view, and the bottom figure shows a cross section through the center of the column, also illustrating the bore. (Reprinted with permission from Motley et al., 2015. Copyright by ASCE.)

In the two-dimensional GeoClaw model, the column was incorporated into the computational domain through the topography term B(x,y) on the right-hand side of Eqs. (2) and (3). Values for B(x,y) are set to a very large constant value, h_c, in the region of the column and to 0 elsewhere. This prevents water from overtopping the area, thus simulating a column. Setting h_c to a very large value also made all four side walls of the square column be more “vertical” in the model since they are represented by steep slopes arising from B=0 (outside the column) to B=h_c (inside the column). The coarsest level grid had a resolution of 0.02 m by 0.02 m and covered most of the computational domain; the finest mesh near the column was 0.25 cm by 0.25 cm, which corresponds to a dimensionless wall distance of $y^{+}=\mathrm{140}$.

First, a case without the column was modeled. Figure 2 shows predictions of water level history, measured at 5.2 m downstream from the gate (i.e., at x=11.1, the center of the column; see Fig. 1 for location of the gauge) by the two numerical models and the experiment. In general, both 2-D and 3-D models accurately predict the arrival time of the bore, which is t=3.2 s.

The OpenFOAM model matches the measurement better than GeoClaw with a sharp (but not vertical) slope at the front, a gradually rising surface to the peak near t=8 s, then a downward slope, followed by interactions with the reflected wave from the back wall that creates the second jump in water level at around t=14 s.

OpenFOAM includes water viscosity, which diffuses sharp discontinuities. In contrast, solving the nonlinear shallow-water equations with an initial discontinuity yields a shock wave (discontinuity) propagating to the right as a vertical bore front followed by a region with constant water depth; as a consequence, GeoClaw slightly overestimates the initial height of the bore front, underestimates the height at t=8 s, and presents the reflected wave as a second sharp discontinuity at t=13.1 s.

Figure 2 Time history of water level at 5.2 m from the gate (center of the column) with the column removed.

At the same location, streamwise (the along-channel direction) components of the velocity at different depths were also predicted. Figure 3 shows time histories of streamwise velocity at nine different distances from the bottom. Note that, since the two-dimensional model is depth-averaged, its predicted velocity is constant with depth. The prediction from the two-dimensional model matches the measurements very well near the water surface, except for the spike at the front, which is better captured by the three-dimensional model. The three-dimensional model underestimating flow velocity near the bottom might be due to our near-wall treatment being imperfect. The velocities in the upper region are also hard to predict well because of air entrained in the water near the free surface as well as the fact that the velocimeter may not be immersed in water at times, causing the measurement to oscillate dramatically.

Figure 3 Time history of streamwise velocity at different distances, d, from the bottom at 5.2 m from the gate (center of the column) with the column removed. Abscissa: time (s). Ordinate: velocity (m s⁻¹).

Figure 4 shows a comparison of total forces on the square column from the experiment, the three-dimensional model, and the two-dimensional model. The force predicted by the three-dimensional model was obtained by integrating the pressure and viscous fluid forces on the surface of the column (see Eq. 17). The three-dimensional model predicts the force very well in terms of magnitude and is able to capture even the small spike near t=4 s. In the two-dimensional model, no hydrodynamic pressure field is available for force prediction. To predict forces from the two-dimensional model, data from the previous case without the column were used instead. The water level, h, and streamwise velocity, u, were first sampled at the center of the footprint of the column that was removed from the domain, to compute the momentum flux, M=hu². As recommended by FEMA P-646 (Applied Technology Council, 2012), the hydrodynamic forces on such a structure can be computed as

$$\begin{array}{}\text{(20)}& F_{\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{d}}\mathit{\rho }\left(h{u}^{\mathrm{2}}\right)b,\end{array}$$

where C_d is the drag coefficient and may be conservatively chosen to be 2.0 as recommended by FEMA P-646 (Applied Technology Council, 2012), F_d is the streamwise component of the fluid forces, ρ is the density of the fluids, h is the water depth, u is the fluid velocity at the location of the structure, and b is the breadth of the structure in the plane normal to the direction of flow. Note that the hu² term in parentheses is the momentum flux, M.

Figure 4 Comparison of measured and predicted horizontal forces on the square column. Sampling frequency is 300 Hz in the experiment and 1000 Hz in both numerical models.

Note that in the experiment or three-dimensional model the water level on the upstream side of the column is different from that on the downstream side of the column. This causes a difference in hydrostatic pressure and thus a hydrostatic force on the column. For this reason, it may be more appropriate to refer to this value as the coefficient of resistance instead of solely as a drag coefficient. Using a drag coefficient of 2.0 overestimates the force by 13 % in general. This is as expected since it is said to be “conservative” according to FEMA P-646 (Applied Technology Council, 2012). Figure 4 also shows that, if a drag coefficient of 1.76 is used instead, the force prediction from the two-dimensional model matches the measurement more closely.

4 The Seaside wave tank model

4.1 The physical experiment

A 1 : 50 scale physical model of part of Seaside, Oregon, adjacent to the Cascadia Subduction Zone (CSZ), was constructed in the Tsunami Wave Basin at the O.H. Hinsdale Wave Research Laboratory at Oregon State University, and a series of experiments were conducted to measure flow velocities and water levels at 31 locations within the model-scale community. For full details on the experiment, one can refer to Park et al. (2013).

The rectangular basin for the experiment is 48.8 m long, 26.5 m wide, and 2.1 m deep. Figure 5 shows the top and side view of the basin. The still water depth at the wave maker is 0.97 m and decreases as it approaches the shoreline. A 0.04 m height (model-scale) seawall was also constructed between all idealized buildings and the shoreline and was parallel to the wave maker. Figures 6 and 7 show the locations of the 31 gauges where water level and flow velocity were measured at a sampling frequency of 50 Hz. The gauges are grouped into four groups – A, B, C, and D (from bottom to top) – and marked by different symbols. Buildings in blue are large commercial buildings like hotels and hospitals. All red buildings are of the same size and represent small commercial buildings. Buildings in yellow are residential structures and are also all the same size.

Figure 5 Top view and side view of the basin.

In the experiment, the piston-type wave maker was designed to generate an initial wave with a wave height of approximately 0.2 m (model scale) at the lower horizontal section of the basin; this is equivalent to 10 m at full scale, which corresponds to a 500-year CSZ tsunami for this region (Tsunami Pilot Study Working Group, 2006). Note that this is not a solitary wave but a long single-peak wave. Experimental measurement of the wave maker speed was fit with a Gaussian function of the form

$$\begin{array}{}\text{(21)}& s\left(t\right)=A{e}^{\mathit{\beta }(t-t_{\mathrm{0}}{)}^{\mathrm{2}}},\end{array}$$

with β=0.25, t₀=14.75, and amplitude A=0.51. Equation (21) was used as input to generate numerical wave in current simulation. The experiment was repeated many times with identical initial conditions. Data from multiple trials were averaged to obtain the results presented here to smooth out stochastic features of the experiment, more details of which were presented in Park et al. (2013).

Figure 6 Layout of all buildings and gauges in the experiment: blue, large hotels or commercial buildings; red, smaller commercial buildings; yellow, residential structures. +: gauge A1–A9; ∗: gauge B1–B9; ⊗: gauge C1–C9; ×: gauge D1–D4.

Figure 7 Four different subsections and layout of gauges.

4.2 Setup of numerical models

4.2.1 OpenFOAM model

In the three-dimensional OpenFOAM model, a numerical wave basin was developed to simulate the experiments. It was built at the model scale instead of full scale to exclude scaling effects. This facilitated the comparison between the numerical model and the physical experiment.

Table 1 OpenFOAM boundary conditions for the current numerical model (for the fixedValue boundary condition, a constant value of 0 is used for the model in this study).

To generate the required waves, a numerical wave generator was previously developed in OpenFOAM (Motley et al., 2014), and it was validated against available data from a pair of experiments. Two steps are taken by the numerical wave generator to simulate the wave-generating procedure of a piston-type wave maker. First, a short subsection of the wave basin adjacent to the wave maker is modeled. This step is conducted with the wave maker as the reference frame, eliminating the need for a moving mesh, and fluid is forced to enter the domain at the wave maker's speed from the end of the domain that is opposite to the wave maker, to simulate the movement of the wave maker. A time-varying acceleration vector field is also embedded in the solver to compensate for the non-inertial frame. The second step is to map all field data in this domain (the generated wave) to a full model of the basin with the mapFields utility in OpenFOAM after the wave maker stops moving. Further simulations can then start from here.

One disadvantage of the three-dimensional model is that it requires heavy computational resources. Even with four dual-eight-core 2 GHz Intel Xeon e5-2650 machines (64 total processors), it was not possible to model the entire basin. Instead, the entire domain was divided into four different subsections of equal width to predict flow parameters at different groups of gauges (see Fig. 7). For clarity, only the onshore domain is shown in the figure; however, the numerical domain spans the entire 48.8 m from the wave maker to the back wall of the basin. For each simulation, approximately 60 million cells were used, and the solver was run in parallel with the 64 processor cores mentioned above for ∼10 days (including wave generation), which is equivalent to a total CPU time of ∼640 days.

The boundary conditions for each boundary in the numerical wave basin are listed in Table 1. The term all walls and floor in the table includes the bottom, side walls, two end walls, and surfaces of internal buildings. Another term, atmosphere, refers to the upper boundary of the computational domain. A zeroGradient boundary condition specifies zero gradient on the boundary. A fixedValue boundary condition sets the value of a quantity to a constant specified value on the boundary. The velocity field on a wall is set to 0. An inletOutlet boundary condition is identical to the zeroGradient boundary condition if the flux is out of domain but is switched to apply a fixedValue boundary condition if the flux is into the domain. The pressureInletOutletVelocity condition at the top of the domain is essentially identical to a zeroGradient boundary condition in our current model. On all walls and floor, p_rgh is defined such that there is zero flux, using the fixedFluxPressure boundary condition, while the atmosphere was defined with a uniform reference pressure p₀ using the totalPressure boundary condition:

$$\begin{array}{}\text{(22)}& p_{rgh}=\left\{\begin{array}{ll}{p}_{\mathrm{0}},& \text{for outflow}\\ p_{\mathrm{0}}-\frac{\mathrm{1}}{\mathrm{2}}|\mathbf{U}{|}^{\mathrm{2}},& \text{for inflow}\end{array}\right..\end{array}$$

Here p_rgh is pressure subtracted by static pressure ρgh, where ρ is the water density, g is the gravitational acceleration, and h is relative depth under the initial free surface. The turbulence quantities near solid walls are obtained with wall functions that model them as functions of distance from the boundary.

Centers of the first layer of cells near the wall are chosen as positions in the log-law region of the boundary layer where the wall functions are applied. A kqRWallFunction boundary condition can be expressed as $\frac{\partial k}{\partial n}=\mathrm{0}$ for k on a wall, where n is a unit normal vector to the wall. An omegaWallFunction boundary condition provides a wall function for the turbulence specific dissipation, ω, with default model coefficients E=9.8, κ=0.41, and C_μ=0.09. It is computed with

$$\begin{array}{}\text{(23)}& \mathit{\omega }=\sqrt{{\mathit{\omega }}_{\mathrm{vis}}^{\mathrm{2}}+{\mathit{\omega }}_{\log }^{\mathrm{2}}},\end{array}$$

where ω_vis is the value of ω in the viscous region, and ω_log is the value of ω in the logarithmic region (Menter and Esch, 2001). The nutUSpaldingWallFunction boundary condition for ν_t is used for smooth walls. It computes a continuous ν_t profile to the wall based on Spalding's law (Spalding, 1961), which is essentially a unified law of the wall which works for the viscous sublayer, buffer layer, and the logarithmic region in a boundary layer.

The initial condition for α_water is set to 1 for cells where there is water at the beginning and to 0 for the rest. The initial values of U and p_rgh were zero since the flow is initially at rest. Although the fluid is at rest at the beginning, a small value of the turbulent kinetic energy k must be “seeded” in the domain, because the production term in the governing equation of the turbulent kinetic energy k is zero and thus will produce no turbulence if initially k is 0.

In the 3-D model, the turbulent kinetic energy k is estimated from the cross-shore velocity u₁, with a factor of 1.25 to take into account the fact that turbulence is three-dimensional (Scott et al., 2005): $k\approx \frac{\mathrm{1.25}}{\mathrm{2}}{u}_{\mathrm{1}}^{\prime \mathrm{2}}$. The velocity fluctuation $u_{\mathrm{1}}^{\prime }$ is computed from $I=\frac{u^{\prime }}{U}$, where I is the turbulence intensity, $u^{\prime }=\sqrt{\frac{\mathrm{1}}{\mathrm{3}}(u_{\mathrm{1}}^{\prime \mathrm{2}}+u_{\mathrm{2}}^{\prime \mathrm{2}}+u_{\mathrm{3}}^{\prime \mathrm{2}})}$, and U can be chosen as wave celerity in this case. This approach is the same as Svendsen (1987) and Lin and Liu (1998). Several choices of initial turbulence intensity were tested, and an turbulence intensity of 1 % is chosen. For the specific dissipation rate, ω, $\mathit{\omega }=\frac{\sqrt{k}}{l}$ is used, where l is the turbulent length scale and is set to 7 % of the hydraulic diameter of the channel-like computational domain, according to Pope (2000).

Two computational meshes with refinement focused on different regions were used during the simulation to mitigate computational demand. The first mesh was used in the first phase, from the beginning of the simulation to the time when the wave almost started to break. In this mesh, all buildings were removed from the domain, leaving only the flat bottom of the wave basin, which reduced the number of cells needed onshore significantly, allowing for use of a much finer mesh offshore. The mesh size was approximately 0.08 m × 0.08 m × 0.01 m (length × width × height) near the wave maker and was gradually reduced to 0.08 m × 0.08 m × 0.004 m onshore due to changes in topography. Note that the mesh cells onshore seem to have large aspect ratio but there is no water onshore at all during this time period. The second mesh was used until the end of simulation. The buildings were added into the domain, and very fine mesh was generated around the the onshore bathymetry. The mesh size was 0.3 m × 0.015 m × 0.035 m near the wave maker and refined to 0.0075 m × 0.0075 m × 0.0025 m near the flat bottom of the onshore segment of the basin and at the edges and corners of the buildings. It should be noted that this was only the size of a structured background mesh, which was further refined by a factor of 2 and deformed by a mesh tool, snappyHexMesh, from OpenFOAM near the buildings and seawalls to make the mesh accurately represent the complex and irregular geometry of the boundaries. Simulation results from the end of the first phase were mapped to the second phase, and the simulation continued. This strategy is similar to the dynamic AMR feature in the 2-D GeoClaw model. Here, however, statically refined meshes were used instead of dynamically refined grids as used in the 2-D GeoClaw model. The average Courant number across the entire computational domain during these simulations is approximately 0.01. While this is considerably low for a typical analysis, this is due to the fact that grid sizes vary by several orders of magnitude.

4.2.2 GeoClaw model

With GeoClaw, it is possible to model the entire basin. Thus, the computational domain is a 48.8 m by 26.5 m rectangle. The geometry of the basin bottom and built environment is described by topography files of different resolution, which specify B(x,y) on the right-hand side of Eqs. (2) and (3). A typical wall time for one simulation is approximately 6 h with a single core in an Intel® Core™ i7-4790 CPU processor, which means the CPU time is also 6 h (0.25 day). Note that the computational resources required by the GeoClaw model are only $\mathrm{0.25}÷\mathrm{640}\approx \frac{\mathrm{1}}{\mathrm{2500}}$ of what is required by the three-dimensional OpenFOAM model in this study.

To generate tsunami waves in GeoClaw, user-defined time-varying boundary conditions can be specified at the inlet of the computational domain, based on data for the wave maker speed s(t) in the physical experiment. The data from the physical experiment can be fit quite well with Eq. (21). However, the way we imposed velocity boundary conditions at a fixed location rather than having a moving boundary, we found better agreement with the observed wave at several offshore wave gauges by setting A=0.6 in Eq. (21), which was therefore used for all simulations.

The AMR feature of GeoClaw was used, with a mesh size for the base-level grid of 0.5 m (corresponding to 25 m in full scale) in both cross-shore direction and along-shore direction. The term cross-shore is used to refer to the direction that the wave propagates from the wave maker to the structures onshore, while the direction perpendicular to the cross-shore direction is referred to as the along-shore direction. The mesh is refined in the nearshore region up to four levels, with specified refinement ratios: 4 from level 1 to 2, 5 from level 2 to 3, and 2 from level 3 to 4. The finest mesh in the domain with this setup for AMR is 0.0125 m by 0.0125 m (corresponding to 0.625 m in full scale) and eventually covers the entire onshore region. The desired Courant number is set to 0.9 to guarantee the stability of the explicit numerical scheme.

One thing to be noted is that, for both numerical models described above, all coastal structures, including different types of buildings and the seawall, are assumed to be undamaged and thus fixed and rigid during the inundation.

4.3 Comparison of flow parameters

The predicted free-surface elevation, cross-shore velocity, and corresponding momentum flux from the two numerical models will be compared and discussed in this section. All experimental data in this study were provided by the NTHMP Mapping and Modeling Benchmarking Workshop: Tsunami Currents (University of Southern California, 2015), and descriptions of the physical experiments to gather the data are provided by Park et al. (2013) and Rueben et al. (2011).

Figure 8 Time histories of surface elevation, cross-shore velocity, and momentum flux at some selected gauges along line A (note that ranges of y axis are different in different subplots).

Figure 9 Time histories of surface elevation, cross-shore velocity, and momentum flux at some selected gauges along line B.

Figure 10 Time histories of surface elevation, cross-shore velocity, and momentum flux at some selected gauges along line C.

Figure 11 Time histories of surface elevation, cross-shore velocity, and momentum flux at some selected gauges in group D.

Gauges were positioned as shown in Figs. 5–7. Ultra-sonic surface wave gauges were used to measure the free surface. The bore front propagation speed was obtained by analysis of imagery gathered by two high-resolution video cameras located above the wave basin (Rueben et al., 2011). Fluid velocity measurements were acquired by acoustic Doppler velocimeter (ADV) only after peaks; air entrainment in the bore at and shortly after the initial impact rendered the ADV measurements inconsistent in repeated trials (Park et al., 2013). Park et al. (2013) then assumed that the propagation speed and fluid velocity at the bore front are equal and fit a second-order polynomial to that value and ensemble-averaged ADV measurements in this region.

Time histories of the free-surface elevation, cross-shore velocity, and corresponding momentum flux at selected gauges are shown in Figs. 8–11. After the peak (initial impact), there appears to be a significant drop in discrepancies between modeled and measured water level and fluid velocity; therefore, the discussion that follows will separately compare the results before and after the peak.

4.3.1 Onshore time series near initial impact

Water level amplitude by OpenFOAM and arrival time by both OpenFOAM and GeoClaw agree fairly well with measurements at many of the gauges in groups A, B, and C, but GeoClaw underestimates the amplitude at many gauges. These differences reflect the challenge of modeling a turbulent and rapidly varying bore front. An additional factor is that the gauges in groups A, B, and C are placed along straight lines, representing roads within the community, whereas those in group D are set behind buildings. As a consequence, flow around group A, B, and C gauges is dominated by flow in the cross-shore direction, while flow around group D gauges is more complex and challenging to model.

Fluid velocity experimental values derived by optical means are significantly lower than the modeled OpenFOAM and GeoClaw velocity in many of the 16 cases presented in Figs. 8–11. This is because the optical measurement of the bore front is not necessarily representative of flow velocity. Here the animation of GeoClaw numerical results was analyzed to obtain estimates of 1.3 m s⁻¹ for peak velocity: Fig. 12 showed modeled velocity distributions in the bore at two consecutive time steps in the GeoClaw simulation at gauge A4, illustrating that the modeled maximum fluid velocity occurs at some point behind the bore front.

Momentum flux modeled by OpenFOAM and GeoClaw do not agree well with experimental estimates, due to the discrepancies in fluid velocity estimates, discussed above. This is critical, since momentum flux is often used to compute the tsunami forces on structure, as discussed in detail in Sect. 5.

In summary, predictions near the initial impact are challenging for both models, but the three-dimensional OpenFOAM model performs better than the two-dimensional GeoClaw model because it models turbulence and the variation of velocity with depth.

Figure 12 Velocity distribution in the bore near gauge A4, from the GeoClaw model.

4.3.2 Onshore time series in post-impact region

Water level agreement among both models and the experimental data is significantly improved after initial impact. Note that some gauges are quite far from the shoreline (for example, gauges A6, B8, and C8), where the inundation depth is very shallow compared to the peak value near the shoreline (less than 20 % of the peak value). Even at these locations, however, both numerical models provide reasonable predictions. It is also of interest that, as noted above, GeoClaw predicts a lower bore front propagation speed than OpenFOAM; as a result, arrival of the OpenFOAM bore front agrees well with experiment, but the GeoClaw bore front is significantly delayed at gauges farther inland, such as B8 and C8 (see Figs. 9 and 10).

Fluid velocity measurements by the ADV are more stable after 30 s, and both OpenFOAM and GeoClaw velocity time series agree much better with the experimental data at gauges in groups A, B, and C. Agreement does degrade significantly in group D, especially in the case of GeoClaw; this is no doubt due to the more complicated fluid flow in the group D environment, behind buildings, compared to the relatively simpler cross-shore flow in the street environments of groups A, B, and C (Fig. 7).

Momentum flux from both numerical models are in better agreement with the measurements at most gauges, since water level and velocity agreements are better than in the t<30 s time period.

Figure 13 compares snapshots of the simulation near line A from the two models at three different times. The three-dimensional model provides substantial detail about the complex flow among buildings, including the strong channeling effect along line A, aligned with the street, and among the buildings on both sides of the street. These channeling effects can alter the forces exerted on both sides of that street, so that any differences between OpenFOAM and GeoClaw in modeling such effects may result in different prediction of forces on the buildings.

5 Force predictions from momentum flux

Some representative buildings along line A were selected for preliminary analysis of fluid forces on the coastal infrastructure, as shown in Fig. 14. Building I is one of the two large structures adjacent to gauge A1. It directly faces the shoreline, with dimensions of 0.29 m by 0.78 m by 0.246 m (length in cross-shore direction by length in 10 along-shore directions by height). Building II has dimensions of 0.39 m by 0.39 m by 0.091 m. Buildings III and IV, representing small houses within the community, are identical but placed in different directions, which have a length, width, and height of 0.17 m, 0.26, and 0.154 m.

Figure 13 Snapshots of the simulation near line A, colored by cross-shore velocity, at three different times (from top to bottom): t=25.9 s, t=27 s, and t=28.1 s. (a, c, e) Geoclaw; (b, d, f) OpenFOAM.

In terms of force measurements, the single-column case presented in Sect. 3 was the only dataset available with experimental measurements of wave impact forces on similar structures. Through validation against that data, it was shown that, provided the water height and fluid velocity are properly modeled, the fluid-induced forces could also be properly predicted. This could be generally extrapolated and applied to the Seaside problem, where the only available measured data included flow parameters (water depth and velocity).

Figure 14 Representative buildings along line A.

Figure 15 shows predicted forces in the cross-shore direction from the two models on selected buildings. Note that these forces are normalized by the width of the western (left) wall of the buildings. Since no pressure field exists in the two-dimensional GeoClaw model, the same approach as was used in Sect. 3 is applied here to compute forces on these selected buildings for the GeoClaw model (C_d chosen as 2.0 as well). In this case, note that not all the buildings are removed to get the momentum flux for a specific building. Instead, only the building at the center of which the momentum flux is to be predicted is removed, with all other constructed environment unchanged. This minimizes the influence of removing that building on the flow overall.

Peak values of forces predicted by the GeoClaw model on all buildings are only approximately half of those predicted by the OpenFOAM model, except for building III. This indicates that this approach for predicting tsunami load can be off by as large as 100 % in such a complex scenario where multiple objects are present, although it is recommended by FEMA P-646 (Applied Technology Council, 2012) as an empirical method when only velocity and surface elevation map is available and is prevalent in tsunami inundation problems.

Figure 15 Predicted forces in cross-shore direction on selected buildings (normalized).

6 Conclusion and extensions

In this paper, two different types of numerical models of tsunami inundation were developed and compared. They were first validated by comparing water level, velocity profile, and forces on a single column impacted by a bore from a dam break. Then the two models were used to predict free-surface elevation, velocity, and momentum flux of a tsunami inundation on a model-scale constructed environment. The predicted flow parameters agree well with experimental measurements in the post-impact region at most gauges. During initial impact, however, the two-dimensional GeoClaw model has difficulty in capturing transient characteristic of the flow. The three-dimensional OpenFOAM model can solve this challenge better, albeit at the expense of many more computational resources being required. This is because the variation in the vertical direction is “eliminated” by the integration in the two-dimensional model, while all three-dimensional characteristics of the flow as well as turbulence are modeled by the three-dimensional model. Several primary conclusions can be drawn from this work:

1. The three-dimensional RANS model can predict flow parameters and forces on structures by modeling a subsection of 1∕3 of the width of the entire basin, while the two-dimensional NSWE model can model the entire basin at one time, with much less computational resources. Both models agree well with experimental measurements at most locations considered after the initial impact. The RANS model, however, can provide more details on the flow, especially near the initial impact region.

2. The fluid dynamics in the bore front are transient and turbulent. Thus near the initial impact, prediction of flow parameters and forces is challenging but also the most critical since the flow parameters and forces have maximum value near this point. The three-dimensional RANS model solves this challenge better than the two-dimensional NSWE model but needs many more computational resources.

3. Using the approach recommended by FEMA P-646 to predict fluid forces on structures from the two-dimensional model works well in the simple case of flow around a column but becomes less reliable in a complex constructed environment. Although choosing a drag coefficient of 2.0 is considered conservative, the 2-D model with this value was still seen to significantly underestimate fluid forces (in some cases giving only half of the prediction from the 3-D model as discussed in Sect. 5) because the 2-D model underestimates peak velocities in this complex flow. For this reason, it is recommended that a 3-D model should be used to determine the tsunami loads on structures when possible, which eliminates the necessity of choosing a large safety factor when only flow velocity is available from the 2-D model as done in Ash (2015).

This research compares different characteristics of a two-dimensional model and a three-dimensional model of tsunami inundation with the constructed environment. Challenges in prediction of flow parameters and forces are revealed, and the capabilities of the two numerical models in solving this type of problem are analyzed. A trade-off needs to be made between the two models due to their different levels of accuracy and required computational resources. The comparisons in the current study can provide a reference when choosing between a two-dimensional model and three-dimensional model to predict required information in tsunami inundation.