Finding an Optimum Grid Size for Numerical Simulations of Dam-break Flow Using Open-Access Digital Elevation Models

This work presents numerical simulations for dam-break flow with open-access digital elevation models (DEMs). The objective is to find an optimum grid size for the simulations in terms of accuracy and computational cost. A hypothetical breaching case of Batutegi dam in Indonesia was selected as a case study. Three DEM types with different spatial resolutions, i.e. MERIT-Hydro (∼90 m), ALOS (∼30 m), and DEMNAS (∼8.1 m), were employed. For each DEM, the numerical simulations were performed using four grid sizes, from the lowest to the finest resolution, i.e. 200 m, 90 m, 50 m, and 20 m. All simulations were carried out on a desktop consisting of 10 cores/20 threads with an in-house code NUFSAW2D, enabling parallelization, where double precision arithmetic was used. It was observed that using the computational grid size finer than 50 m produces no longer significant differences of the maximum inundation boundary, maximum depth, maximum velocity, and average depth. This study would be useful for engineers and decision makers, when performing dam-break simulations using DEMs in a framework of an early-warning system.


Introduction
Extreme flood events in the last decade have demonstrated the large damage to people life and economy all over the world. However, assessing flood hazard impacts, one of which by numerical flood simulations, remains a challenging task up to now. Therefore, comprehensive research in flood modeling needs to be continued and updated in order to deepen the knowledge as well as to improve the capability of predicting the flood hazard.Two ofextreme flood events in Indonesia were: (1) the Gintung dam-break case in 2009, where approximately 2 million m 3 of water were suddenly released, 100 people were dead, approximately 100 people were missing, and a total of 10 ha of the downstream area was inundated [1], and (2) the Way Ela dam-break case in 2013, which released approximately 20 million m 3 of water causing one people dead, one people missing, and 32 people injured [2].Albeit smaller, Gintung dam was more hazardous than Way Ela dam because no emergency action plan was comprehensively considered yet. This shows an importance of a dam-break risk assessment that may be conducted through numerical simulations. The state-of-the-art in carrying out a dam-breaksimulation is to usedigital elevation model (DEM). Over the last decade, the availability of high-resolution DEMs has increased through remote sensing and field surveys [3], making accurate predictions of flood modeling more feasible. While highresolution DEMs (finer than 30 m) are usually not free of charge and in the absence of direct topography measurement, open-access DEMs, e.g. Multi-Error-Removed Improved-Terrain (MERIT) Hydro, Advanced Land Observing Satellite (ALOS), and DEMNAS may be an option for flood modeling. MERIT Hydro is a ~90 m resolution DEM that has been developed utilizing global hydrography datasets and multiple inland water maps, thus containing flow direction, flow accumulation, hydrologically adjusted elevations, and river channel width [4]. ALOS is a ~30 m resolution DEM developed utilizing the Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) [5]. DEMNAS is an Indonesian DEM product with a resolution of ~8.1 m, which is a combination of various sources, i.e. IFSAR (~5m resolution), TERRASAR-X (~5m resolution), and ALOS PALSAR (~11.25m resolution) with Mass point stereo-plotting data [6]. Essentially, DEM resolution plays a significant role in the accuracy of flood inundation modeling. The computational grid size for a numerical flood model can be set equal to or even finer than the DEM resolution in order to achieve accurate results. Nevertheless, the computational cost may become tremendously expensive. Moreover, the finer grid sizes may not always exhibit significant increases in accuracy. While yet the effects of DEM resolution on the flood inundation results are still rarely investigated, it is therefore of uttermost importance to study and understand the trade-off between the computational grid size and the computational cost.
In this work, MERIT Hydro, ALOS, and DEMNAS are used for the flood simulations of a hypothetical breaching case of Batutegi dam in Lampung Province, Indonesia. The focus of this study is on investigating the effect of various grid sizes on the accuracy and computational cost of the numerical dam-break simulations for different DEM resolutions. For each of theaforementioned DEMs, the simulations are carried out using different grid sizes, i.e. 200 m, 90 m, 50 m, and 20 m, to observe the differences of flow characteristics including depth and velocity. Note that one should use ALOS and DEMNAS carefully for hydraulic simulations because they are of a digital surface model (DSM), whereas MERIT Hydro is of a DEM, to which corrections with several hydrologic and hydraulic features, i.e. river lines, have been applied, thus being appropriate for flood modeling. In our experience (unpublished), the uses of DSM for some cases of hydraulic simulations have ledto incorrect results. However, this issue is not discussed in detail in this work as no anomaly was observed for the numericalresults using ALOS or DEMNAS compared with the ones using MERIT Hydro.
This study uses NUFSAW2D (Numerical simUlation of Free surface ShAllow Water 2D), an inhouse code developed by the first author of this paper. NUFSAW2D solves essentially the 2D shallow water equations with some extensions to turbulence and 2D/3D non-hydrostatic models. NUFSAW2D is currently supported with hybrid OpenMP-MPI (Open Multi-Processing/Message Passing Interface) parallelization, while the GPU (Graphics Processing Unit) version is under development. NUFSAW2D has been proven accurate and efficient, see the documentations [7][8][9][10][11][12][13][14][15][16]. Note that for the sake of brevity, the mathematical and numerical formulations of NUFSAW2D are not given here, thus the interested readers may refer to the aforementioned publications. The final objective of this paper is to find an optimum grid size for dam-break simulations using open-access DEMs. Here, optimum is in the sense that using grid sizes finer than the limit size will no longer give significant differences of the flow characteristics, i.e. inundation boundary, depth, and velocity.

Description of Case Study
The study area, Batutegi dam, located in Lampung Province, Indonesia, with a coordinate of 5°15'17.32"S and 104°46'48.86"E. At the downstream of Batutegi dam, there is Way Sekampung dam constructed in 2016, look at figure 1. As no cascade dam-break case is considered here, the domain for the Batutegi dam-break simulations is limited to the upstream part of Way Sekampung dam. In other words, once the water from Batutegi dam arrives at Way Sekampung dam, the numerical simulations stop. Several data are prepared for this study such as DEMs, bed roughness map, and hydrograph value. For the sake of simplicity, only DEMs and hydrograph data arebriefly discussed in the next two sections.

Digital elevation models
Three DEMs obtained are shown in figure 2. The domains are focused on the area between Batutegi and Way Sekampung dams, being approximately 17x12 km. The long sections along the river line as well as three cross-sections (at the upstream, middle, and downstream parts of the river) are presented in figures 3 and 4, in which the bed elevations have been adjusted using the lowest bed value among the three DEMs at the upstream part of the river (Station 0 m). Figure 3 shows that the bed elevations of ALOS and DEMNAS are relatively higher than those of MERIT Hydro. Figure 4 shows that at the upstream part of the river, the channel cross-section of MERIT Hydro is wider than the others, while both ALOS and DEMNAS show similar cross-sections. At the middle and downstream parts of the river, all DEMs exhibit similar cross-sections but MERIT Hydro has lower bed values than the others. This may be due to the DSM type of ALOS dan DEMNAS, where as MERIT Hydro has experienced the corrections and adjustments of hydrology and hydraulic features.

Technical data of Batutegi dam and breach hydrograph
The technical data of Batutegi dam are briefly presented in table 1. According to [17], Batutegi dam is safe against overtopping during the probable maximum flood (PMF) event. However, the piping failure mechanism is of consideration. Due to the difficulty in computing the piping process even with the help of numerical tools, the piping breach hydrograph in [17] was estimated using the empirical formula of [18], which only gave a final breach shape of a trapezoidal form with a side slope of 1:0.7 (V:H) and an average width of 234 m. However, the mechanism of the breach propagation is not exactly known, and thus, it must be assumed to define the hydrograph value. Yet, it is supposed to follow a linear propagation process, i.e. from the initial phase (a rectangular hole as an orifice flow) to the final phase (a trapezoidal breach as a weir flow), following [19]. An empirical criterion that governs the mechanism change from the orifice to weir flow is employed following [20]. The detail of breach mechanism is not discussed here. The process of the breach propagation toward its final shape is sketched in figure 5.  Figure 6 shows the total outflow hydrograph during the breach event that represents a total value of the spillway outflow (during the PMF event) and the breach outflow (both orifice and weir flow mechanisms). One can observe that the spillway of Batutegi dam gives a peak outflow of 781 m 3 /s, approximately at time = 23 hours, during the PMF event. The water depth above the spillway is approximately 3.83 m (+277.83 m). At this point, the initial piping hole is assumed to occur on the dam body at the elevation of +229.5 m. In the next 1.52 hours or at time = 24.52 hours, the final breach shape is reached. This condition causes a peak total outflow discharge of 186,585 m 3 /s, emptying the reservoir in about 2.6 hours from the first hole formation.  Figure 6. PMF inflow and outflow hydrographs (left) and the hydrograph of the total outflow during the breach event (a total value of the spillway and breach outflows) (right)

Comparison of Accuracy
First, it is important to note that this paper does not aim to present the comparison of accuracy between MERIT Hydro, ALOS, and DEMNAS because no measurement data are available yet. Instead, the focus on investigating the effect of the computational grid size on the numerical results. Therefore, the comparison of accuracy using different computational grid sizes is presented only for each DEM. One can see in figure 7 that using the grid size of 200 m for the ~90 m MERIT Hydro significantly produces larger inundation areas than the others. This is clearly shown at the downstream location (at the upstream part of Way Sekampung dam), where the inundation areas are significantly larger, for example, than the ones computed with the grid sizes of 50 m and 20 m. The water depth distributions using 200 m grid size especially at the centerline of the river are also different and tend to be smaller than the others. Refining the computational grid size from 200 m to 90 m has produced a remarkably different result. This can be seen at the downstream location, where the inundation areas become significantly smaller. Refining the computational grid size from 90 m to 50 m still exhibits remarkable differences at the downstream location. However, they are not as significant as the differences produced when refining the grid size from 200 m to 90 m. Interestingly, it is observed that using 50 m grid size seems not to give significant differences in inundation areas with the ones computed using 20 m grid size. This shows a remarkable decrease of the inundation boundary, in which the use of 200 m grid size forevery DEM resolution seems to overestimate the flood effect. The grid size refinements from 90 m to 50 m and 50 m to 20 m exhibit similar behaviors for each DEM, being 0.87-0.89x and 0.95x, respectively. This clearly shows that a further grid size refinement of less than 50 m will not produce any significant improvement of the maximum inundation boundary.
With respect to the maximum depth, the refinement of the grid size from 200 m to 90 m for every DEM brings the results to an increase of the depth value about 1.10-1.13x. Reducing the grid size from 90 m to 50 m gives a maximum increase of 1.06x, while doing so from 50 m to 20 m only results in a maximumincrease of 1.04x. This indicates that refining the grid size less than 50 m does not produce any significant difference for the maximum depth value. In terms of the maximum velocity, the grid size refinement from 200 m to 90 m has tremendously increased the value, being 1.49x, 1.48x, 1.44x, for MERIT Hydro, ALOS, and DEMNAS, respectively. However, with the further grid size refinements, the differences in the maximum velocity value become smaller. For example, the refinement from 90 m to 50 m produces the differences in the range of 1.10-1.15x, while the refinement from 50 m to 20 m is in the range of 1.07-1.13x. Again, one may see that the grid size refinement of less than 50 m does not yield any significant improvement for the maximum velocity value.
Finally, the accuracy comparison was presented in terms of the average depth. An interesting fact is observed here, where the simulations with MERIT Hydro have experienced the smallest increase of the average depth of 1.22x by the grid size refinement from 200 m to 90 m. Meanwhile, for both ALOS and DEMNAS, the same grid refinement has increased the average depth value 1.40x and 1.31x, respectively. However, the refinements from 90 m to 50 m and 50 m to 20 m bring each DEM to experiencing similar increases of average depth value in the range of 1.11-1.13x and 1.04-1.05x, respectively. This again shows that with 50 m grid refinement, one may not expect any significant improvement for the average depth value.

Comparison of Computational Cost
In this section, a comparison of computational cost is presented by taking as a base value the average CPU time required for the simulations using 200 m grid size. All simulations were performed using 10-cores/20 threads. Note that the average CPU time is quite similar for all DEMs with the same grid size. Therefore, the comparison shown in table 3 has represented the average CPU time for MERIT Hydro, ALOS, and DEMNAS. Another important parameter that has been taken into consideration is the Courant-Friedrichs-Lewy (CFL) condition as it governs the stability of the explicit numerical computations used in this study. In order to present a fairly-representative comparison (for the computational efficiency with respect to the total number of grids), only one constant time step value was selected for all numerical simulations, namely 0.5 second, which is the value that still maintains the CFL number of less than 1 for the simulations with 20 m grid size.  In the perspective where one constant time step value of 0.5 second is used for all simulations, the increase of the total number of grids causes a tremendous increase of CPU time, showing a nonlinear relationship. One of the reasons is due to the parallelization overhead cost. In another perspective, where one may use a larger value of time step (especially for the simulations with grid sizes larger than 20 m), refining the grid size may be viewed as a worse case, albeit with the help of parallel computation, as it will tremendously increase the actual CPU running time. Therefore, modelers should take this issue into careful consideration, when attempting for a considerably shorterreal running time.

Conclusion
The numerical simulations of the Batutegi dam-break case has been presented using an in-house code NUFSAW2D for MERIT Hydro (~90 m), ALOS (~30 m), and DEMNAS (~8.1 m). It was observed that using the computational grid size of 200 m for MERIT Hydro, ALOS, and DEMNAS seems to overestimate the flood effect by larger inundation areas but may underestimate it by lower average depth and lower maximum velocity. Refining the grid size from 200 m to 90 m for each DEM producedsmaller inundation areas as well as larger maximum depth, maximum velocity, and average depth. A similar behavior was also noticed for the grid size refinement from 90 m to 50 m. Hereinafter, refining the grid size from 50 m to 20 m exhibited no longer significant improvements with respect to the accuracy.
In terms of computational cost, the increase of the total number of grids caused a tremendous nonlinear increase of CPU time. From all the results, it is reasonable to conclude that 50 m grid size is quite optimum in terms of accuracy with relatively lower computational cost for practical engineering purposes in dam-break simulations, when using (open-access) DEMs, not only with low-but also high-resolution maps. As only the rapidly varying unsteady flow is considered in this work, the investigation of the grid size effects for the numerical simulations with the gradually varying unsteady flow would be an interesting topic for future study.