Investigating the Capability of HEC-RAS Model for Tsunami Simulation

This study highlights the simulation of tsunami cases using the freeware HEC-RAS 6.1. The primary aim is to evaluate the capability of the software in performing tsunami simulation due to its standalone computational framework (pre-processing, execution, and post-processing stages), making the modeling process interactive especially for practical engineering purposes. The model accuracy was tested against some benchmark cases of wave propagation, including analytical solutions, laboratory experiments, and field measurements. The results were mainly compared based on the amplitude fluctuations of water surface elevation and velocity magnitude. For the analytical case, the elevation and velocity were accurately computed. Furthermore, sufficiently accurate results were obtained for the laboratory case, where the maximum elevation was properly computed without significant discrepancies. For the field cases, the wave arrival time and the fluctuations of water surface and velocity were also appropriately calculated. The errors produced were assessed using the Root Mean Square Error (RMSE) and the Pearson Product Moment Correlation (PPMC) parameters. The RMSE values between the numerical results and the analytical/observed data were relatively low below 30% and the PPMC values were ranging from 52–99%. Based on these results, it can be stated that HEC-RAS was capable of modeling tsunami propagation. In addition to its eminence, a drawback was found regarding the graphical user interface (GUI) of HEC-RAS for the input of boundary conditions, which is indicated by two issues. First, it was the flow hydrograph (instead of the velocity value) used for the boundary condition, which affected the numerical computation for the momentum terms of the shallow water equations. Secondly, the limitation of the time interval was considered to cause inaccuracy result to a certain extent. These findings will be beneficial for the coastal engineering community and the continuous development of HEC-RAS.


INTRODUCTION
A tsunami can be triggered by several phenomena, such as earthquakes, landslides, and even volcanic eruptions, and can travel trans-oceanic distances at high speeds without losing a significant amount of energy (International Institute for Geo-Information Science and Earth Observation (ITC), 2005).Tsunami is different from other waves, such as tidal and wind-generated waves, in terms of wavelength and period.Typically, its wavelength and period are distinct, ranging between 100 -500 km with a period of up to 1 hour.In contrast, tidal wave has a wavelength of up to 15 km, while windgenerated wave is only 100 -200 m, each with a period of up to 12 hours and approximately 10 s, respectively.Tsunamis in the open ocean can be classified as shallow water waves, where the wavelength is much longer than the water depth.During their propagation process, the velocity com-ponent is integrated with the water depth, causing the waves to slow down as they approach the shore.Consequently, their height rapidly grows, reaching up to 3 m in 90 seconds (Dias and Dutykh, 2007).The accurate prediction of tsunami height remains a challenge as it is difficult to estimate the source parameters (Sannikova et al., 2021).
The propagation of tsunami can be studied and understood in several ways, namely historical tsunami events, physical modeling, and numerical modeling.Among these methods, numerical modeling has become a powerful tool with the rapid growth of computational resources over the past few decades.In general, the model can express some important physical parameters of the flow propagation phenomena, and the output accuracy depends on how well the phenomena can be integrated into the numerical schemes (Thomas and Dwarakish, 2015).Predicting tsunami inundation becomes even more challenging since the onshore tsunami behavior is directly related to the beach slope and the structure's shape (Prasetyo et al., 2019).The numerical models mostly used for tsunami simulations can be classified into 2D and 3D models (Hajihassanpour et al., 2019).Although the 3D models are expected to produce more accurate results compared to 2D, they are quite impractical for large-scale cases, as the computational process could take considerable time.Therefore, the 2D models are still preferable in the scope of tsunami simulations.
This study summarizes the advances in tsunami simulations based on some works from the 1960s -2020.The numerical tsunami simulation was initiated in Japan and developed by Aida (1969) and Aida (n.d.) to simulate the 1964 Niigata and 1968 Tokachi-Oki tsunamis.Goto and Ogawa (1982) also presented numerical simulations for tsunami problems using the finite difference method (the Leapfrog scheme).Imamura (1989) further developed the model, which solved the shallow water equations using the bottom friction term.Shuto et al. (1990) also extended the model to NAMI-DANCE, by considering the tsunami source either from a defined ruptured waveform or water surface fluctuation, as an input.Liu et al. (1995) developed COMCOT for 1D/2D tsunami modeling using the finite difference method, while Numerical modeling of 3-D long wave runup using VTCS-3, au-thor=Titov, VV and Synolakis, CE (1996) numerically simulated long wave runup using the VTCS-3 model.Extensive developments of tsunami models were made in the 2000s.One notable contribution was Tsunami-HySEA, introduced by Castro et al. (2005).This model utilized a highorder finite volume scheme and was employed for many tsunami prediction cases with complex bathymetry and overland flow.Furthermore, Geo-Claw was developed at the University of Washington and pioneered by George (2006) to simulate tsunami propagation using a high-resolution shock-capturing finite volume method.
The MOST model was developed by Wei et al. (2008) for tsunami modeling covering their generation, propagation, and inundation onto dry land.Kim et al. (2009) developed pCOULWAVE to solve the weakly-dispersive, turbulent Boussinesq-type equations using a high-order finite volume model.In contrast, Roeber and Cheung (2012) designed BOSZ to specifically solve the near-field tsunami and extreme swell waves (hurricanes) problems.Tolkova (2014), used an open-source version of the MOST model to simulate a tsunami by applying a splitting method that reduces the 2D propagation problems to two 1D problems.In addition, an open source code (FUNWAVE-TVD) based on the Boussinesq modeling was used to compute the coastal inundation for the wave propagation (Shi et al., 2012).A Fortran-based source code (TUNAMI-N2) was used by Oishi et al. (2015) to forecast the near-field tsunami inundation for the 2011 Tohoku-Oki earthquake case.With an open source code (sam(oa) 2 ), the simulation of the tsunami propagation resulting from the 2011 Tohoku earthquake in Japan, was carried out using billions of computational grids (Meister et al., 2016).In a recent study, the tsunami propagation was modeled using an in-house code (NUF-SAW2D) to solve the 2D non-hydrostatic shallow water equations through an artificial viscosity technique, where all the results were consistent with the observed data (Ginting and Ginting, 2020).
Most of the aforementioned models are in the form of source codes that may require external compilers for execution and external post-processor tools to visualize the results.Another common way for practical engineering purposes in modeling tsunami is the use of software that is already supported with a Graphical User Interface (GUI).One such is MIKE-21, which is a commercial tool developed by Danish Hydraulic Institute (DHI).With this tool, the users can interactively work from the pre to post-processing stages by using a well-developed GUI, but with a required subscription fee.As a developing country, Indonesia has limited funds for practitioners to purchase commercial tools for their projects.Therefore, freeware supported with GUI is preferable.Another option for tsunami modeling is DELFT3D (in 2D mode).This is a freeware that utilizes several (external) open-access tools with GUI support, such as DELFT DASHBOARD, MUPPET, and QUICKPLOT for pre and post-processing stages (Deltares, 2022).
This study evaluated the capability of HEC-RAS 6.1 for tsunami simulations, with HEC-RAS partic-ularly chosen for two reasons.First, it is a standalone freeware supported with GUI from the preprocessing (building meshes), execution (fullydynamic shallow water solver), to post-processing stages (visualization), therefore making the modeling process interactive for the users.HEC-RAS has been frequently used for hydrodynamic simulations and some of the recent works are noted below.It was used to simulate the flood inundation on a complex and highly anthropogenically altered topography (Shustikova et al., 2019).It was also used to simulate the hydrodynamic rainfallrunoff processes at a basin scale, highlighting its potential for the water engineering community in the innovative study field (Costabile et al., 2020).Moreover, HEC-RAS has been employed to investigate the effect of land cover on the tsunami overland flow propagation (Adityawan et al., 2021).
The second reason for choosing this model is its constant process of update and continuous improvements, both in terms of advanced numerical scheme and technology for the computational speed-up.The evaluation of the capability of HEC-RAS, specifcally for tsunami modeling, is of benefit to practitioners.Despite its eminence, there was a limitation in the model's GUI for the input of boundary conditions.Therefore, these findings can contribute to the continuous development of HEC-RAS.

HEC-RAS 6.1
The governing equations used in HEC-RAS 6.1 are based on two laws, namely mass and momentum conservation.In terms of mass conservation, the flow is assumed to be incompressible.Regarding the momentum equation, the vertical velocity is assumed to be negligible as the horizontal scale is much longer than the vertical one.Therefore, the vertical velocity is neglected, and the pressure is assumed to be hydrostatic.The differential forms of both equations are as follows: where h is the water depth, ∇ is the differential operator, V denotes the velocity vector in x and y directions (u and v), t represents time, q is the source/sink flux term, zs is the water surface elevation, fc is the Coriolis parameter, k is a unit vector in the vertical direction, vt is the eddy viscosity tensor, τb and τs are the bottom shear stress (computed using the Manning formula) and the wind surface stress vector, respectively, while R is the hydraulic radius.It is worth noting that this study used a solver for the shallow water equations with an Eulerian-Lagrangian approach for the advection term (SWE-ELM), where the Coriolis term, turbulent mixing model, and wind surface stress are neglected.
HEC-RAS 6.1 uses a combination of the finite difference and finite volume methods on unstructured meshes to solve governing equations in spatial discretion.When the grid is orthogonal, the finite difference method is used to calculate the normal derivative, while for non-locally orthogonal grids, the normal derivative is split as the sum of a finite difference and finite volume approximation.The mesh can be calculated using sub-grid bathymetry.The main idea of this approach is to use high-resolution bathymetry data to calculate the sub-grid levels on coarse grids, allowing large time steps (Sehili et al., 2014).For the temporal discretization, the Crank-Nicolson method can be used to weigh the variables' contribution at time step n and n+1.The wetting and drying process of multiple cells can be implemented in a single time step by using a semi-implicit scheme.The other flow properties, such as hydraulic radius, volume, and cross-sectional area are pre-computed from the bathymetry and saved in the computational grid cells.This approach is suitable for various applications.The details of the numerical scheme of HEC-RAS are presented by Brunner (2021).
To generate a proper computational mesh, each cell is required to capture the bathymetry and the changes in water surface slope and velocity.Therefore, the selection of a proper grid size is quite tricky, as it can significantly affect the accuracy of results and computational cost.Moreover, the grid size also influences the choice of the time step value, as a smaller grid size often demands a smaller time step value.Alternatively, the time step in HEC-RAS 6.1 can be determined by setting the maximum and minimum Courant numbers.
Subsequently, the model adjusts the proper time step automatically.The Courant number equation is defined as follows: where C is the Courant number, V stands for the velocity, ∆T and ∆X are the time step and the average cell size, respectively.
The HEC-RAS uses the GUI to define the computational domain and establish the boundary condition, with automatic mesh generation.In this study, the rectangular domain was simply set for each test case according to its domain size.A single mesh size was subsequently specified in the GUI and HEC-RAS to automatically create the computational meshes, thereby producing rectangular, uniform meshes.However, when a nonrectangular (free-form) domain is used, HEC-RAS will generate uniform meshes, except for some areas near the domain boundaries, where the other shapes are automatically determined by HEC-RAS.
The flow boundary is established by setting the boundary line either inside or outside of the domain.There are three types of boundary lines, namely water surface elevation, discharge value, and normal depth.To run HEC-RAS, it is necessary to set at least one flow boundary line in the domain.Once a flow boundary line is established, the other boundaries are automatically set as walls, preventing any flux from entering or leaving the domain.In contrast, some tsunami models require at least two boundary lines to define an inflowoutflow process, with the outflow boundary line usually treated with an absorbing condition.As for the initial condition on the computational domain, three types are available, namely single water surface, dry condition, and restart file.In cases where the second option is applied, the cells remain dry until the water flows into the domain.Meanwhile, the last option is used to establish values for an entire simulation by providing (different) water surface elevation for every cell in the model.Therefore, the other flow parameters can be computed based on the sub-grid bathymetry approach.For all test cases in this study, the initial water level was automatically defined from the initial value of the stage boundary.The bed elevation is automatically treated as a dry bed in HEC-RAS when it is higher than the initial value.However, it is worth noting that the specification of velocity as a boundary/initial condition is not currently available in HEC-RAS, as detailed in Sections 3 and 5.

Focus and Objective
In general, this study aims to investigate the capability of HEC-RAS 6.1 in modeling tsunami cases.
To demonstrate its eminence, the numerical results were compared based on three benchmark categories, namely analytical solution, laboratory experiment, and field measurement.Foremost, a sinusoidal wave case was selected as a representative of the analytical case.Subsequently, a case of laboratory benchmark experiment, namely a tsunami runup onto a complex 3D beach, was presented.This is one of the benchmark problems in The Third International Workshop on Long-Wave Runup Models held at Wrigley Marine Science Center, California (ISCE, 2004).
Furthermore, the 2011 Tohoku tsunami event was selected to represent the field measurement category.This event was recorded in two different locations, namely Hilo Bay, Hawaii and Tauranga Harbor, New Zealand, and their results were compared.These cases were the benchmark problems presented in the National Tsunami Hazard Mitigation Program (NTHMP) workshop held in Portland, Oregon (at University of Southern California , USC).The subsequent section explains the accuracy of the numerical results, which were quantified and analyzed with statistical parameters.
In addition to its eminence for tsunami modeling, the limitations of HEC-RAS was studied with respect to its GUI for the input of boundary condition.First, this limitation relates to the type of boundary condition.Currently, the model does not support the input of velocity value as a boundary condition, despite that in some tsunami cases, their usage is mandatory in conjunction with water elevation.Although this problem can be mitigated using flow hydrograph value from the GUI, it might decrease the accuracy of the numerical results to a certain level.The second issue that may significantly cause inaccuracy is the time interval limit for the boundary condition set by the GUI, which cannot exceed 1 s.For some certain tsunami simulations, a smaller time interval is required to achieve better accuracy.This paper explores this limitation by examining a laboratory case of a solitary wave on a conical island conducted by Briggs et al. (1995).

RM SE
To assess the error produced by HEC-RAS 6.1, two statistical parameters were used.The first parameter is the Root Mean Square Error (RMSE), which illustrates the spread of residuals between the regression data of the numerical results and the observed data.The bigger the RMSE value, the larger the error between the two sets of data.The RMSE formula is expressed in Eq. 4.
where N is the total number of data considered, Dcomp is the numerical result, and Dobs denotes the observed data.Secondly, the Pearson Product Moment Correlation PPMC was used to indicate how well the numerical results and the observed data are related.In other words, the PPMC value represents the linear relationship between two sets of data or measures the strength of the relationship between two variables.The PPMC is expressed as: (5) There are three typical values of PPMC.A value of 1 indicates that for every positive increase in the numerical results, there is a positive increase of a fixed proportion in the observed data.Conversely, a value of -1 means that for every positive increase in the numerical results, there is a negative decrease of a fixed proportion in the observed data.
A value of 0 indicates that there is no positive or negative increase, or in other words, the two sets of data are not related.

Case 1: Idealized Inlet with Sinusoidal Wave
The first test case in this study was based on the thesis of Ginting (2011).The objective was to simulate a sinusoidal wave propagation in an idealized inlet with a frictionless, rectangular flat channel, measuring 5 km in length and 1 km in width.The Manning coefficient was set to 0.00001 s m -1/3 , which is close to zero since HEC-RAS cannot operate with a zero manning coefficient.The properties of the wave were set to an amplitude of 2.5 m and a period of 3 hours and 12 hours.These properties were applied as a boundary condition along the 1 km boundary line, while the others were set as wall boundary conditions.The water was initially at rest with a depth of 10 m and assumed to be at an elevation of +0 m.The analytical solutions for the wave surface and velocity are as follows: where a, L, x, t, η, and u, are the amplitude, channel length, distance, time, wave surface, and wave velocity, respectively.Note that v is zero.
The domain was discretized with a grid size of 20 m.The total simulation time was set to 72 hours.However, the comparison between the analytical and numerical results is presented in this study as the simulation was considered stable after 24 hours.The computation time was around 11 minutes, running with Intel® i7 CPU @1.30GHz with an 8.00 GB RAM Laptop (these specifications were also used for the other test cases).For this case, a comparison was made between the numerical and analytical results at the center of the domain.Fig. 1 shows that HEC-RAS 6.1 can produce accurate results.The comparison of water level showed no discrepancy for both simulations.However, there was an insignificant discrepancy in velocity, specifically at t = 12 hours.In other to address this, there was an approximate 11-minute time shift to capture the maximum velocity.As shown in Fig. 4, the numerical results obtained at all stations were consistent with the benchmark data.The first incoming wave was captured at 14.4, 14.6, and 14.9 seconds at gauges Ch-5, Ch-7, and Ch-9, respectively.The maximum elevation was also properly simulated with merelyinsignificant discrepancies, although the numeri- The first 10 s could not be captured appropriately due to the existence of initial water disturbances in the wave tank (Zhang and Baptista, 2008).Similar discrepancies in water level during the first 10 s were also observed by Nicolsky et al. (2011) and Hou et al. (2015).Even the 3D model used by Kakinuma (2008) could not compute them appropriately.
Fig. 5 shows the visualization of the wave propagation between 0 -17 seconds.It can be observed that the wave remained relatively stagnant for the first 10 seconds and started rising at 15 s, inundating the beach area between 15.5 s and 16 s, where the maximum elevation was reached at around 17 s.Subsequently, the water level gradually decreased until all reference points were dry again.During this period, the complex wet-dry mechanism was noticeable, while the stability and accuracy of HEC-RAS 6.1 in simulating it were demonstrated.every 60 seconds from 7 hours to 13 hours (6-hour window).The computation process for this case was completed in approximately 1.5 hours.
Fig. 8 illustrates that HEC-RAS 6.1 effectively simulates the water level at Tide Station.The first wave arrival was accurately captured at around 8.2 hours, while the fluctuation pattern had a positive correlation with the benchmark data.The water surface elevation was also well predicted at 8.4 and 8.6 hours, respectively.This indicates that the lowest water surface elevation can be appropriately estimated.However, the numerical result for the maximum elevation was overestimated by 0.32 m compared to the observed data of +1.14 m.Overall, HEC-RAS 6.1 produced similar results to the observed data.This can be seen from the complex water surface fluctuations from 7 -13 hours that can be properly simulated with merely non-significant discrepancies.
For the comparison of velocity, the observed data were collected at a 6-minute interval, while the numerical results were plotted with a 1-minute interval.The first wave arrival can be accurately predicted in HAI1225 station at approximately 8.2 hours.According to Fig. 9 (a), the numerical model appropriately predicted the first wave arrival, with approximately 0.25 m s -1 at 8.2 hours.However, the subsequent results tend to be overestimated, with 0.3 m s -1 observed at 8.3 hours.The amplitudes still show similar patterns compared to the observed data.The gap between the numerical results and observed data tends to increase over time.Nevertheless, the values after 10.5 hours had a better prediction.erally in line with the observed data at HAI1226 station.This is indicated by the maximum velocity amplitude that can be properly predicted at approximately 1 m s -1 .However, there are still differences between the numerical results and the observed data.Three phenomena can be noted in this regard.Firstly, the numerical results detected the first incoming wave by capturing it 3 minutes faster than the observed data.Nonetheless, this discrepancy is still tolerable in practical applications.Secondly, the numerical results showed a In this context, the numerical results at 8.3 hours only satisfied 50% of the observed velocity.Lastly, there were discrepancies in predicting the velocity fluctuations from 7 -10 hours, although they improved after 10.5 hours.These discrepancies may be attributed to the dispersion phenomena originating from the ocean or tidal current that was not numerically modeled.Moreover, discrepan-cies in the water surface and velocity results were also found in the similar previous work of Arcos and LeVeque (2015).
In addition to the comparison across the three stations, it is also important to visualize the velocity propagation during the simulation time to understand the distribution over the domain.First, the results for 8.2 hours were presented when the first incoming wave at Tide Station was detected.According to Fig. 10 (a), there was an increase in velocity with approximately 1 m s -1 at the toe of the breakwater.A similar increase was also observed at Tide Station, with a notable surge of 1 m s -1 near the coastline (river estuary).The breakwater remained dry at this point, with no overtopping water.
The next visualization was presented at 8.4 hours, when the minimum water surface elevation was observed, as shown in Fig. 10 (b).It can be seen that the water above the breakwater is still dry.
The maximum velocity at the toe of the breakwater increased up to 3.2 m s -1 , and the high-velocity distribution started spreading towards the outer part of the breakwater.Meanwhile, the velocity at Tide Station and near the river estuary now decreased toward zero, with water subsequently overtopping the breakwater.This is illustrated in Fig. 10 (c), when the maximum water surface elevation was observed at 8.6 hours.Currently, the breakwater is inundated by water, with the velocity above the breakwater ranging between 1.4 -2.1 m s -1 .Meanwhile, the maximum velocity at the toe of the breakwater significantly increased, reaching approximately 3 m s -1 .
The velocity at Tide Station and near the river estuary increased to 1.2 m s -1 and 1.8 m s -1 , respectively.The breakwater was dry at 10.5 hours and the velocity at the toe of the breakwater decreased between 0.8 -1 m s -1 .Meanwhile, an increase in velocity was observed in Tide Station, at approximately 1.8 -2.2 m s -1 .Therefore, the velocity was concentrated at the outer part of the breakwater, ranging between 0.8 -1 m s -1 .It is evident that HEC-RAS 6.1 can stably simulate the wet-dry mechanisms, even for small depth values of 0.01 -0.001 m.The maximum velocity distribution for the entire simulation time is visualized in Fig. 11.For clarity purposes, the scale used ranges between 0 -5 m s -1 .The maximum velocity is mainly concentrated in front of and at the toe of the breakwater, with values ranging between 2.5 -3.5 m s -1 , which significantly increased at the top of the breakwater.The velocity distribution tends to be more uniform at the outer part of the breakwater than at the inner.This distribution pattern is in accordance with the numerical results presented by Ginting and Ginting (2020) and Arcos and LeVeque (2015).m for velocity comparison.The boundary condition or incident wave was given along the northern shoreline part.Similar to Case 3, the incident wave values in Case 4 indicated only the tsunami signal as the effect of the tide had been removed from the raw data (at University of Southern California , USC).The time window was set from 10 to 20 hours, with 10 seconds time step.The Manning coefficient was set to 0.025 s m -1/3 , which was similar to the previous case.The domain was sketched in Fig. 13, while it took approximately 1 hour to successfully run this case.
According to Fig. 14, the numerical model appropriately captured the first wave arrival at all stations for both water level and velocity.The first wave arrival was captured at approximately 13.3 hours.Meanwhile, amplitude fluctuations showed a similar pattern, where the lowest and highest elevations can be properly simulated at each station.
Although the results were generally in accordance with the benchmark data, there were still some discrepancies.Firstly, the numerical results tend to overestimate the observed data.For instance, at 17.8 hours the numerical result overestimates the water level by 0.15 m in Tug Berth station and the velocity by 0.4 m s -1 at ADCP point.The water level was also overestimated by 0.18 m at 18.5 hours in the Moturiki station.
Secondly, a time delay was detected after 16 hours, which is most noticeable in Tug Berth, Sulphur Point, and ADCP stations.Although a delay of up to 7 minutes was detected in Sulphur point around 17.1 hours, these discrepancies are still considered tolerable for practical engineering purposes.ranging from 0 -0.7 m s -1 , while it remained relatively uniform at the outer region.Fig. 15 (c) shows the velocity distribution when the water elevation is minimum.Contrary to the former, the velocity distribution tends to be more uniform inside the harbor, with approximately 0.4 -0.8 m s -1 at the harbor entrance.
Finally, Fig. 15 (d) shows the presentation of the visualization when the velocity is maximum.The maximum velocity is mainly concentrated at the harbor entrance alley, with a maximum value of approximately 1 ms-1.Interestingly, HEC-RAS 6.1 is again capable of simulating the wet-dry mechanisms properly with the small values of depth ranging between 0.01 -0.001 m, similar to Case 3.This reveals that the model is sufficiently accurate for real tsunami simulations with complex bathymetry and wet-dry mechanisms.

Drawback on the GUI for Boundary Condition
Having presented the capability of HEC-RAS for tsunami modeling in the previous section, the limitation of this model in relation to its GUI is subsequently discussed.As previously mentioned, a solitary wave on a conical island case was selected (Case 5).According to Fig. 16, the basin size was 25 x 30 m, and the conical island was placed at the center of the domain.The Manning coefficient was set to 0.00001 sm -1/3 .For the purpose of comparing the experimental data and the numerical results, three gauges were selected, namely P-3 (6.82; 13.05) m, P-6 (9.36; 13.80) m, and P-9 (10.36; 13.80) m.The domain was discretized into a 250 x 300 grid or 75,095 cells.The computational time was set to 0.1 seconds, with a total simulation time of 20 seconds.The left side boundary was set as the flow boundary, while the other three sides were treated as walls.The incident wave was specified as a solitary wave, with the water elevation and velocities respectively defined as: where Ae is the amplitude,He denotes the constant water depth, and Te is the time the crest enters the domain.In this case, these parameters were set to 0.032 m, 0.32 m, and 2.45 seconds, respectively.It is worth noting that while specifying the boundary condition, the GUI of HEC-RAS 6.1 only supports the input for the water surface elevation (stage hydrograph) and not the velocity.Therefore, the flow hydrograph (which is a multiplication of the water depth, channel width, and velocity at the left side boundary) was used to represent Equations ( 8) and (9).The computational process took approximately 21 minutes to completely run the model, and the incident wave propagation is presented in Fig. 17.At t = 9 seconds, the run-up arrived first, producing a wet and dry mechanism on the island.Also, the maximum water surface elevation was reached during this period.Meanwhile, at t = 11 seconds, the incident wave started colliding in front of the island and subsequently refracted because the water depth changed and started to propagate around the island.Therefore, after the collision of the wave at t = 13 seconds, the second run-up was generated, and the wave subsequently propagated around the island to the entire domain.
According to Fig. 18, the numerical results showed discrepancies at all observed gauges.The numerical model captured the first wave arrival time at approximately 0.6 -3 seconds faster, with a longer period than the experimental data.In terms of the maximum elevation, the numerical results at Gauges P-3, P-6, and P-9 were sufficiently in line with the experimental data, where the measured discrepancies were only around 0.002 m.The wave elevation decreased over time as the maximum elevation was exceeded.
After thoroughly observing the GUI of HEC-RAS 6.1, it is important to emphasize two points.First, the discrepancies shown in Fig. 18 may be attributed to the use of flow hydrograph as the boundary condition.Since the HEC-RAS performs computation from cell center to cell center, similar to most 2D finite difference/finite volume models, the use of flow hydrograph instead of water level and velocity along the boundaries may affect the momentum calculations at the cell faces, and subsequently the other parameters within the internal domain.When water level and velocity are used at boundaries, the water depth and velocity vector for the subsequent time step will directly be advanced in the numerical scheme, exactly following the given/known boundary condition.However, this is not the case when the flow hydrograph is employed, as there is another procedure to convert such a hydrograph back to depth and velocity before it can be advanced for the subsequent time step.Nevertheless, it is likely that this problem is not the sole reason for the discrepancies observed in Fig. 18.
The second problem worth addressing is the limitation of the time input interval when specifying the flow boundary condition.The GUI limits the interval time to 1 second, whereas the implementation of the flow boundary condition (Equations ( 8) and ( 9)) for this case requires a level of precision of 0.01 seconds to properly capture the maximum water surface elevation.Therefore, the accuracy of the flow boundary values is only 1/100 of what it should be.
This could be specified as the main cause of the discrepancies between the numerical results and observed data in Fig. 18.
This condition is illustrated in Fig. 19, where the comparison of elevation and flow data produced a 1-0.01-second level of precisions.The maximum elevation was 0.352 m when a 0.01-second time interval was applied, while it was only 0.340 m when 1-second was used.
A similar condition applies to the flow hydrograph data, where the discrepancy for the maximum discharge was up to 0.67 m 3 s -1 .
To support the hypothesis, several trials were conducted by slightly adjusting the shapes of the flow hydrograph.For example, (a) the value of Te was adjusted between 2 -3 seconds and (b) Ae value was adjusted between 0.033 -0.034 m.For the first attempt, the numerical results could only capture the first arrival time accurately, but the maximum water elevation was underestimated.Conversely, for the second attempt, the numerical results were able to predict the maximum water elevation, but not the first arrival time.Therefore, it can be con-cluded that the discrepancy in this case was not caused by the core solver of HEC-RAS 6.1, but by the limitation of the time interval provided by its GUI.

Accuracy Assessment
In this section, Table 1 presents the errors of the numerical results in the forms of RMSE and PPMC.
For Case 1 at T = 3 hours, the respective values of RMSE and PPMC for the water elevation were 0.08 m and 0.99, while 0.05 m s -1 and 0.91 were obtained for the velocity.The maximum error magnitude for this case ranges between 5 -8%, which can be considered accurate.The results indicate a strong relationship between the analytical and numerical results, with a correlation coefficient of up to 99%.When T = 12 hours is applied, the RMSE and PPMC values for the water elevation were 0.01 m and 0.99, respectively, while 0.01 m s -1 and 0.97 were obtained for velocity.The average error magnitude produced was only 1%, which is lower than in the previous case.In addition, the correlation coefficient ranged between 97 -99%.
In Case 2, the numerical results showed a good agreement with the benchmark data, with all sta- possibly tend to exist between 9.5 -10.5 hours, hence the maximum and minimum magnitudes for this range cannot be properly captured by HEC-RAS 6.1.Based on the PPMC value of 0.57, the correlation coefficient for the water surface elevation between both data was estimated to be 57%.The comparison of velocity at HAI1225 station showed that the values of RMSE and PPMC were 0.17 m s -1 and 0.68, respectively.The former indicates that the average error magnitude of the velocity was only 17%, which is lower than that of water surface elevation.A better PPMC value of 0.68 is shown at this station, indicating that about 68% of both numerical results and observed data tend to move in the same direction.Therefore, the correlation coefficient between both data was 68%.Meanwhile, the values of RMSE and PPMC at HAI1226 station were 0.23 m s -1 and 0.52, respectively.The first value is slightly higher, while the second is slightly lower than the ones at HAI1225 station.
In Tauranga Case, the RMSE parameter values at Tug Berth, Sulphur Point, Moturiki, and ADCP Point were 0.08 m, 0.10 m, 0.11 m, and 0.15 m s -1 , respectively.The lowest error was recorded at Tug Berth station, with only 8% of the average error magnitude.This value indicates a high accuracy between the observed data and the numerical results.Meanwhile, the respective values for the PPMC parameter were 0.78, 0.68, 0.79, and 0.80.Therefore, the relationship between the observed data and the numerical results can be considered strong, with a correlation coefficient of up to 80%.
Finally, the RMSE values for the Conical Island were much smaller compared to the other cases because the GUI requires a higher level of precisions for the input of the boundary conditions.The most accurate results can be seen in P-3 with RMSE and PPMC reaching 0.0047 m and 0.86, respectively.Although the results are considered sufficiently accurate, the discrepancies in Fig. 18 are still significant.

CONCLUSION
This study highlights the results of tsunami modeling using the freeware HEC-RAS 6.1.Based on the findings, the following conclusions can be drawn: HEC-RAS 6.1 is capable of performing the tsunami simulation and successfully producing stable numerical values.The results were generally considered accurate, as demonstrated by the relatively low errors below 30% and the strong average relationship of 52 -99% between the computed and observed data.This can be seen in Case 1 -Case 4, where the amplitude fluctuations generally showed a similar pattern, while the benchmark data, the maximum water level, and velocity were appropriately predicted.
The discrepancies between the numerical results in Case 5 can be attributed to the limitation of the HEC-RAS GUI for the input of boundary condition, which can be classified into two problems, namely the boundary condition type and the time interval limit.Although the use of a flow hydrograph (instead of water level and velocity) as boundary condition can decrease the accuracy to a certain level, it was hypothesized that the discrepancies in Case 5 were more likely caused by the time interval limit when specifying the flow boundary condition.Therefore, the GUI limits the interval time to 1 second, whereas the implementation of the flow boundary condition for Case 5 requires a level of precision of 0.01 seconds.Addressing this limitation could be beneficial for the future development of HEC-RAS.
HEC-RAS 6.1 is a standalone freeware supported with GUI from the pre-processing (building meshes), execution (fully-dynamic shallow water solver), and post-processing stages (visualization).Therefore, it can make the modeling process interactive for the users.Considering that the aforementioned limitations will be fixed in the future, HEC-RAS can be regarded as a reliable tsunami modeling tool for the coastal engineering community in the scope of practical engineering purposes.

Figure 1
Figure 1 Case 1: Comparison of wave surface and wave velocity between analytical and numerical results at (a) t = 3 hours and (b) t = 12 hours

Figure 3
Figure 3 Case 2: Incident wave for the boundary condition

Figure 4
Figure 4 Case 2: Comparison between the numerical results and experimental data

Fig. 9
Fig. 9 (b)  shows that the numerical results are gen-

Figure 6 Figure 7 Figure 8
Figure 6 Case 3: Bathymetry, domain area, and reference points

Figure 9
Figure 9 Case 4: Velocity comparison between numerical results and observation data: (a) HAI1225 and (b) HAI1226

Figure 16
Figure 16 Case 5: Conical Island domain plan

Table 1 .
Comparison of error values