Adaptive Mesh Refinement for Dam-Break Models using the Shallow Water Equations

The 2D shallow water equations are a common tool for the simulation of free surface fluid dynamics in civil engineering. However, the nonlinear structures of the equations’ straightforward implementations lead to numerical problems, such as spurious oscillations and unphysical diffusion. Therefore, this research compared several strategies to overcome these problems, using various finite element formulations and combinations of stabilization methods and mesh options. The accuracy and performance of numerous approaches are examined on models of dam-break in one and two space dimensions. The analytical solution checks the numerical, derived shock wave heights and velocities for the 1D classical benchmark. The result showed that streamlined diffusion and shock capturing stabilization deal with the classical problems of spurious oscillations and numerical diffusion but still indicate similar problems locally in the vicinity of steep fronts and shock waves when used on fixed meshes. As adaptive meshing is the most promising method to deal with such situations, several concerned options are examined in detail. It is important to fine-tune the method to the model’s needs, i.e. to adapt the maximum number of mesh refinements, the indicator functions, and the starting mesh. The use of adaptive meshing techniques leads to accurate solutions for the usual parameter range in 1D and 2D, requiring less computational resources than simulations on fixed meshes. Meanwhile, meshing reduces the model size of the 2D dam-break model adaptive by almost one order of magnitude and the execution time by a factor of 20.


INTRODUCTION
The shallow water equations (SWE), also known as Saint-Venant equations, are fundamental for modeling fluid flow in shallow open fluid systems. The SWE has been applied for a wide field of applications, such as open channels (Chaudhry et al., 2008), floods (Bermúdez et al., 1991), coastal currents (Brocchini and Dodd, 2008), tsunamis (Liu et al., 2009;Takase et al., 2011) and estuaries (Li and O'Donnell, 1997).
The major reason for its application is that the dimensionality is reduced from three to two or one. In the real world, hydraulics could be described by the 3D Navier-Stokes equations. For most applications, which include turbulence, these equations have to be extended to deal with the complexity of fluid flow. Even for simple geometrical set-ups, the general equations require large computational resources, storage, and execution time. Therefore, from the viewpoint of the modeler, the reduction of the dimensionality is a crucial step. The SWE can be written as follows: with total water depth, water height above reference height, velocity vector, and acceleration due to gravity denoted by H, η, u, and g (Takase et al., 2011). The equations are derived from the volume and momentum conservation principles formulated on depth-averaged variables. The derivation is based on several assumptions, such as fluid incompressibility, hydrostatic pressure distribution in the vertical direction, depth-averaged values for all properties and variables, small bottom slopes, and no density effects from variable fluid density or viscosity. Furthermore, the eddy viscosity is much larger than the molecular, and the atmospheric pressure gradient can be ignored.
Generally, the solution of the nonlinear SWEs (1) and (2) may suffer from severe instabilities because straightforward modeling without stabilization, using finite differences or finite element techniques, leads to spurious oscillations. This led to the proposal of various stabilization schemes, such as the introduction of an artificial viscosity, which appears in an additional term on the left side of equation (Chen et al., 2013).
where artificial viscosity is denoted by v. This equation is determined in analogy with stabilization methods for the advection-diffusion process, where similar numerical problems arise simulating almost sharp chemical and thermal fronts. The artificial viscosity method prevents oscillatory errors but introduces a numerical diffusion, which is not presented in the original equations and real systems. The originally proposed method was inconsistent with a grid-dependent artificial viscosity modified using finite volumes by Ginting (2017) to obtain second-order accuracy.
There is a vast literature on consistent stabilization methods for finite elements, differences, and volumes approaches, which all have advantages, disadvantages, and limitations. The streamline diffusion method (Donea and Huerta, 2003), also termed 'Streamline Upwind Petrov Galerkin' (SUPD) (Hughes, 1979), is popular in general finite element discretization. It modifies the diffusivity and effectively smoothens the solution while increasing errors near the sharp fronts.
Therefore, combining streamlined diffusion with spurious oscillations is convenient for correcting errors near the front using diminishing layers (SOLD) methods, such as shock-capturing diffusion (Augustin et al., 2011). According to John and Knobloch (2007), the combination of SUPG and SOLD is consistent. The following formulas describe the weak formulation for the 1D case.
H and u denote the test functions for water height and velocity. δ stream serves as a switch to include or dismiss streamlined diffusion and τ SU P G is a tuning parameter. The formulation for two velocity components for the 2D case has to be specified analogously.
However, there are still problems with the transition of sharp fronts independent of the stabilization method. This is because the general step function cannot be represented by solutions that are only defined in a discrete space. When the front is located between mesh points, the numerical solution cannot match the analytical solution.
Such an error cannot be reduced by a better representation of physical processes but only by mesh refinement.
Solutions are therefore needed with inconsistent and consistent stabilization by analyzing the finite element functions of a different order. This study demonstrates that local or global mesh refinement effectively provides solutions to practical problems that are too large to handle. Adaptive meshing is the method that yields the highest accuracy in relation to the use of computer resources.

METHODS
The various approaches are examined for 1D and 2D dam-break scenarios, which are convenient test problems for the SWE. The use of dam-break models for benchmarking software for the solution of the SWE equations dates back to ancient times (Biscarini et al., 2010;Erpicum et al., 2010;Baghlani, 2011;Peng, 2012;Duran, 2015;Ginting, 2017;Buttinger-Kreuzhuber et al., 2019;Vichiantong et al., 2019;Putri et al., 2020). A rough examination of the literature reveals that several set-ups of dam-break models have been used for benchmarking with differences in dimension in the conceptual model and values.
The dam-break models used for benchmarking are based on a highly simplified concept that excludes the appearance of the breaking wall in a simple setting. The simulation starts at time t=0 with the steep step of the water table at the interface between the reservoir and the backwater region. Based on the highly idealized step situation, the SWE simulates the development of the water table for t>0. The simple set-up includes two water heights, including one higher value for the reservoir (h 1 ) and a lower (h 0 ) for the backwater. In a 1D model, the interface between the reservoir and backwater is at the position x 0 , where the dam and the initial water table jump are located as shown in Figure 1.
The existence of an analytical solution, first derived by Stoker in 1957 using the characteristics method, makes the dam-break problem attractive. It consists mainly of two waves of a different kind, namely a simple wave moving into the reservoir and a shock wave in the backwater region. The sketch in Figure 2 shows the characteristic shape of the water table at a time instant.
The water depth of the shock wave h 2 can be calculated as a function of the initial heights h 1 and h 0 , as shown in Equation 4 (Wu et al., 1999).
Note that equation differs from the formula given by Wu et al. (1999) but delivers the same solutions if X and Y are taken as X = h 0 /h 1 and Y = h 2 /h 1 . The nonlinear equations are solved numerically to compute the celerity of the shock wave in a second step. The water table in the simple wave region is given by the quadratic function in equation: In 2D, the dam-break models used for benchmarking differ in the geometry of the model regions and the initial conditions. A quadratic model region with a 200 m side length and 75 m long partial break of a straight dam was used in several benchmark studies (Biscarini et al., 2010;Baghlani, 2011;Vosoughifar et al., 2013;Jalalpour and Tabandeh, 2014). Pilotti et al. (2010) focused on the appearance of the reservoir waves after a partial break of a straight dam in several other geometries. Erpicum et al. (2010) constructed Lshaped model regions while conducting a dambreak study.
The radial dam-break model is simpler due to an unrealistic circular dam, which surrounds the origin of the coordinate system. The mathematical analysis can reduce this problem to 1D using a cylindrical coordinate system. Despite the unrealistic features, a convenient test case for the numerical models in 2D was used to quantify errors induced by meshes in Cartesian coordinates (Erpicum et al., 2010;Baghlani, 2011;Remacle et al., 2006;Pilotti et al., 2010;Jalalpour and Tabandeh, 2014).
In the sequel, one 1D and two 2D models were examined. The interval x ϵ [0, 1], and x 0 = 0.5 were chosen as the region and position of the dam in the 1D model with a physical unit and model region extension of 1 m. The models can be conceived as valid for other physical dimensions if all length values are non-dimensionalized by their length as a reference value. Furthermore, when not mentioned otherwise, the reference models in 1D and 2D have a flood height and backwater height of 1(h 1 = 2, h 0 = 1) at a time interval of 0.1 s. As the set-up becomes more challenging if the backwater height is small, the errors in a parameter study are explored with a modified ratio h 0 /h 1 . In the 2D set-up, the unit square is chosen as the model re- Figure 3. Comparison of numerical and analytical results for shock wave heights and velocities for (1) fixed initial step height and varying backwater depth and for (2) fixed backwater depth and varying initial step height gion, and locate the interface at the circle around the origin of radius r = 0.5. This makes the dam geometry a quarter of a circle.
The 1D reference model in the numerical solution has 100 linear elements and a regular mesh size of 0.01. In the refined 1D model, 400 elements and a mesh size of 0.0025 were used, while in the 2D model, it was 0.02, which correspondents to 2500 linear rectangular elements (mapped mesh). The refined 2D model has a mesh size of 0.01, corresponding with 10000 rectangular elements, while the double refined mesh consists of 40000 elements. Furthermore, adaptive meshing was utilized for triangular meshes.
These simulations deal with a time interval in which the waves do not reach the boundaries because when it does, additional numerical errors are expected. The additional errors depend on the type of condition, which is dissimilar to the wall boundary conditions used in all model runs in this study. All simulations were run us COM-SOL Multiphysics (2021), a software operated using pre-defined physics modes and by specifying the partial differential equations in a finite element formulation (pde mode). The software is applied in many application fields, mostly in engineering. The model for the SWE was set-up using the weak formulation in the pde mode. Only the most recent version of COMSOL contains a predefined mode for the shallow water equations that were not used in this study.

1D Dam-Break
The dam-break model was simulated with various values for initial step height h 1 and backwater depth h 0 at reference values of 1 m, using the refined mesh model (400 elements) with consistent stabilization. Figure 3 shows the shock wave height and velocity results in dependence of h 0 and h 1 . Both plots show two graphs for constant h 1 and h 0 , which changed between 0.05 and 1 m, as noted on the x-axes.
Both plots depict the numerical results compared with the analytical solutions, with the latter obtained using Equations 5 and 6. These variations are not problems for the numerical solver to produce the correct values of the analytical solutions. A closer look shows deviances for low ratios h 0 /h 1 .
The model was run with lower ratios to further explore the numerical method's limitations. Figure  4 depicts the relative errors for depth ratios down to 0.01 and shows that the error increases strongly when the ratio h 0 /h 1 is lower than 0.05. At that threshold, the relative error lies at approximately 1% and increases to more than 5% and 8% for the shock wave velocity and height when the ratio decreases to 0.01. Low backwater depths represent a challenge for the numerical approach.
In the sequel, the front development is presented in plots of water height at different time instances. Figure 5 depicts simulated water tables obtained using different numerical approaches, demonstrating some previously mentioned known issues. The upper plot illustrates the problem of oscillations appearing without applying stabilization to produce an unacceptable solution. The middle plot compares the results of model runs using consistent and inconsistent stabilization. The grey graphs depict solutions obtained using an artificial viscosity, while the colored graphs are determined by combining streamlined diffusion and shock-capturing stabilization. The inconsistent stabilization results in a numerical diffusion that smoothens all gradients. This numerical method makes sharp slopes less steep, which is also unacceptable.
The bottom plot of Figure 5 illustrates the effect of using different element orders, whereby the grey and colored graphs result from using linear and quadratic elements. The figure indicates that the quadratic elements better resolve the steep front because the spatial resolution of the numerical approach is finer than those of linear elements. However, a more detailed analysis reveals that at the front quadratic elements are more prone to oscillations and instabilities than their linear counterpart. Therefore, linear elements are preferable because a better suppression of numerical diffusion is paid for by a higher number of degrees of freedom (DOFs). DOFs of the refined model, at execution times of 9 and 19 s, respectively.
A comparison of the results on the two equidistant meshes, (a) and (b), shows that the finer mesh suf- Figure 6. Height output for the 1D dam-break model; detailed view near the shock wave front at times t=0.09 and 0.1, detailed view of shock wave front: (a=grey) reference mesh, (b=black) refined mesh, (c=blue) adaptive meshing, (d=red) fine-tuned adaptive meshing fers less from numerical diffusion, while the sharp front is better resolved. However, mesh refinement seems to have almost no effect concerning over-and undershooting errors in the backwater and along the shock wave. This is because the amplitude of the oscillations is almost the same.
Additional runs were performed using the adaptive mesh refinement process to explore the role of mesh refinement in more detail. The results were reported using default settings (c) and run with adjusted, fine-tuned parameters (d). The element growth rate was increased from 1.7 to 3, with a rise in the maximum number of mesh refinements from 2 to 6. The adaptive mesh runs had DOFs mean values of 346 and 705, with execution times of 2:12 and 7:29 (min:s).
The spatial resolution plays the most important role, while mesh refinement reduces the numerical diffusion to reproduce the sharp front movement. Results obtained with adaptive meshes show the least numerical diffusion, suppressing oscillations with an increase in element growth rate and maximum mesh refinements. However, for better solutions, the number of DOFs needs to be increased for higher storage and execution time requirements. The costs and benefits of adaptive mesh refinement are explored in detail for the 2D dambreak benchmark.

2D Dam-Break
This section thoroughly examines the 2D dambreak models with a quarter circular dam. Figure  7a shows the initial state and the wave propagation following the break event on the right. The This corresponds with the simple wave region of the 1D case with the penetration of a shock wave in the opposite direction into the backwater region. In contrast to the 1D analytical solution, the shock is not a pure step function. The transition from the simple wave to the shock wave region appears quite abrupt.
The 2D solutions along the two diagonals of the square model region were analyzed to determine the error along the main diagonal of the longitudinal wave movement. However, the results are unacceptable when the model is run without stabilizing spurious oscillations in the 1D benchmark. Figure 8a demonstrates that numerical diffusion induced by inconsistency is more pronounced in the 2D case than in the 1D simulation. Here the results with the artificial viscosity approach are compared to the streamlined diffusion and shockcapturing stabilization. The consistent stabilization technique in Figure  8b shows the effect of mesh refinements on the front development in the longitudinal direction. The reference mesh was refined twice using a grid spacing of 0.01 and 0.005 m to produce steepening fronts, which are effective on the left side at the transition between simple wave and shock regime. The higher accuracy, the more expensive the price of computer resources. The reference mesh has 7803 DOF, and the simulation took 43s execution time, which indicates that the refined model needs 30603 DOF and 10:23 min execution time. The double refined mesh has 121203 DOF and needed 1 h 23:41 min for execution. This indicates that the execution time increases by 14 due to the first mesh refinement and 8 in the second refinement.
The differences between the solutions and the errors become more obvious when plotted on a cut line along the transverse diagonal of the model region. The results obtained are compared to determine the inconsistency and consistency of the stabilization, as shown in Figure 9a.
The middle plot shows that the effect of mesh refinement is better than in Figure 8. The output of the first refinement (gray) does not make much difference, while the second steepened the fronts substantially. A big jump in accuracy can be observed between the first and second mesh refinement, which comes at extra costs, as outlined in Figure 9. Several model runs aimed at reducing the computational costs using adaptive meshing were also documented. The elements are modified locally from the initial coarse mesh depending on the current solution. In unsteady models, mesh refinement is also combined with a decrease in time-step. Several options control adaptive meshing, such as the size of the initial coarse mesh. There are different ways to refine an element geometrically, including the provision of an element growth rate to control the coarsening of the mesh. Also, maximum refinements are required to prevent the model from becoming too large and reduce the execution time.
The criterion of which elements have to be refined is based on a local error estimate, whose functions are determined by the user. It is usually based on the change of the dependent variables, such as H, u, and v. Although the variable can be used, it is better to utilize derivatives as error indicator functions (Chellamuthu and Ida, 1995;Grätsch and Bathe, 2005). Two indicator functions are used based on the height H, as shown in 7 and 8: The other is based on the velocities u and v and thus reads (Garcia and Popiolek, 2014): The user also specifies the number of sample points on which the error is estimated, as well as the options for the combination of different variables for the local error estimation. Figure 10 shows the typical sequence of meshes taken during a run with adaptive meshes. Only one segment with a refined mesh can be identified when the two waves are close in the initial phase. Therefore, with the movement of both fronts, two segments arise, separated by a zone with coarser elements.
The bottom plot of Figure 9 compares the results using the double refined mesh with those obtained by the adaptive type, which was run with an initial spacing, maximum refinement, sample points, and error indicators of 0.02, 3 and based on Equations 8. There are no visible deviances on the front flanks, but the center contains an adaptive mesh solution, which is better than the fixed mesh result, showing some unphysical oscillations. Table 1 provides the list of model runs with adaptive meshing, which lists the parameters of the grid adjustment in the first four columns. The 6 th column represents the mean model size of 10 meshes produced during the simulation in addition to a measure of the model size of degrees of freedom (DOFs). The execution time on a common computer with 2.2 GHz processor is given in the 7 th column, while the last measures the error. The calculation error due to the unavailability of the analytical solution for the 2D problem was based on the use of a numerical solution. Therefore, the L2-Norm was used to calculate the difference between the numerical solutions in question, which is the one for the double refined mesh. It turns out that a) the number of sampling points is of minor importance: for coarser meshes, there are no influences (4 5), while for finer meshes (runs 14 15), a higher number of starting points increases the accuracy and reduces the run time. b) the error indicator should be based on the  Figure 10. Development of adaptive mesh, example derivatives, to increase the accuracy of the results by comparing run 3 with 4, and 12 with 15. c) the best error indicator is based on the derivatives of the velocity components d) increase of the maximum number of refinements by comparing runs 15 and 16, which may lead to a significant improvement in accuracy. However, such an improvement is not guaranteed, as indicated by comparing runs 5 and 6. e) it is important to have an initial mesh of appropriate refinement, as shown by comparing runs 5 and 15, as well as 6 and 16.
The improvement due to the tuning exercise, doc-umented in Table 1, indicates that the coarse mesh results in the upper part of the list consist of a mesh, which was increased by a factor of 4. Furthermore, the execution time increases by the roughly same factor, although the gain in accuracy was by more than one order of magnitude.
The comparison of the finest mesh adjustment with double refined mesh solution shows that mesh refinement reduces the DOF by almost a factor of 8. In terms of execution time, the effect of adaptive meshing is more striking because the variable mesh run finished in only a 20 th of the time of the fixed mesh run.

DISCUSSION AND CONCLUSSION
Numerical simulations of the SWEs are likely to become inaccurate when the solution shows steep fronts and shock waves. The behavior of the model implementations under these circumstances is easily studied using 1D and 2D models of dam break.
The results confirmed that straightforward implementations and inconsistent stabilization lead to disturbed solutions and are unacceptable. Several stabilization methods, like streamlined diffusion and shock capturing, can be used to deal with the problems of oscillations and numerical diffusion.
However, implementations with fixed mesh are unable to generally cope with the positioning of sharp fronts or shock waves. This is because the advancement of a shock wave cannot be resolved in a numerical approach based on a fixed mesh. The height at the front will turn into a gradient, which can be conceived as an effect of numerical diffusion accompanied by spurious oscillations at the front position. For such situations, the accuracy can only be improved by refined meshing.
Global mesh refinement is involved high costs on computer resources, such as storage and execution time, which means it has to be performed locally. The use of 2D makes it easy for the computer storage and execution time to be achieved. The process of moving fronts mesh refinement can only be obtained through adaptive meshing, which is a convenient technique because it remains within appropriate limits.
Adaptive meshing was examined using a 2D dam break model, and the simulations demonstrate the importance of fine-tuning the method to suit the model's needs. The importance of turning is that it determines the maximum number of mesh refinements, the indicator functions, and the starting mesh. It also reduced the adaptive meshing model size and execution time by 8 and 20 compared to the fixed globally refined mesh.
Adaptive meshing is an economical and efficient technique used to obtain accurate solutions to problems containing shock waves. The dam-break applications showed an increase in errors with small backwater heights using the numerical approach, making it difficult to deal with wetting and drying sub-regions. Relative to the vast literature on stabilization methods, adaptive meshing as a topic does not seem to have attracted the attention it deserves.
[This page is intentionally left blank]