## Abstract

A wavelet-based local mesh refinement (wLMR) strategy is designed to generate multiresolution and unstructured triangular meshes from real digital elevation model (DEM) data for efficient hydrological simulations at the catchment scale. The wLMR strategy is studied considering slope- and curvature-based refinement criteria to analyze DEM inputs: the slope-based criterion uses bed elevation data as input to the wLMR strategy, whereas the curvature-based criterion feeds the bed slope data into it. The performance of the wLMR meshes generated by these two criteria is compared for hydrological simulations; first, using three analytical tests with the systematic variation in topography types and then by reproducing laboratory- and real-scale case studies. The bed elevation on the wLMR meshes and their simulation results are compared relative to those achieved on the finest uniform mesh. Analytical tests show that the slope- and curvature-based criteria are equally effective with the wLMR strategy, and that it is easier to decide which criterion to take in relation to the (regular) shape of the topography. For the realistic case studies: (i) slope analysis provides a better metric to assess the correlation of a wLMR mesh to the fine uniform mesh and (ii) both criteria predict outlet hydrographs with a close predictive accuracy to that on the uniform mesh, but the curvature-based criterion is found to slightly better capture the channeling patterns of real DEM data.

## INTRODUCTION

Spatially distributed rainfall–runoff modeling can rely on the numerical solution of the shallow water equations – for example, Simons *et al.* (2014); Coon *et al.* (2016); Caviedes-Voullième *et al.* (2020) – and is becoming popular in hydrological modeling (Fatichi *et al.* 2016). Numerical solvers for the shallow water equations are currently robust, efficient, and parallelized, enabling to simulate catchment-scale problems at high resolutions (Xia *et al.* 2019). However, such simulations still entail prohibitive computational costs on a uniform mesh that significantly grows with an increasingly finer resolution (Xia *et al.* 2019). As topography data are increasingly available at higher resolutions than 10 m, it is desirable to find an alternative strategy to reduce computational cost for such high-resolution simulations at the catchment scale (Ferraro *et al.* 2020).

For spatially distributed rainfall–runoff, hydrodynamic and hydrological responses are mainly controlled by small topographic features which can impact runoff generation (Thompson *et al.* 2010) and drive runoff connectivity across different portions of a large domain (Bracken *et al.* 2013). This property of the topography in relation to hydrological processes is referred to as the terrain connectivity or structural connectivity (Bracken *et al.* 2013). A foreseeable strategy is to use a multiresolution mesh, where the small features of the topography are captured at a locally refined resolution on the mesh while efficiently employing coarser resolution elsewhere to reduce operational cost. This can be achieved by an adaptive mesh refinement (AMR) strategy (Berger & Colella 1989), many of which are already established for shallow water modeling (Shashidharan *et al.* 2018, and cited references within). AMR can be applied dynamically in time to adapt the resolution of the mesh at every time step or can be static in time to generate an initial multiresolution mesh. As the flow hydrodynamics induced by rainfall–runoff is primarily controlled by static topographic features, a static AMR strategy can be sought, and is referred to here as *local mesh refinement* (LMR) (Caviedes-Voullième *et al.* 2012).

Applying an LMR strategy to digital elevation model (DEM) data to generate a reliable multiresolution mesh is not a trivial task: on one side, there is the issue of the solution's sensitivity to mesh design; and, on the other, there is no literature evidence that topography-based a priori error estimators exist for hydrological flow processes. Thus far, LMR strategies require ad hoc criteria for deciding local coarsening and refinement of mesh resolution. These criteria require inputs based on either slope or curvature data of the DEM (Caviedes-Voullième *et al.* 2012; Hou *et al.* 2018; Ferraro *et al.* 2020). Slope-based and curvature-based refinement criteria may lead to refinement at different locations; therefore, a comparative study for identifying the implications of such choices for LMR on multiresolution mesh decision is still required. A related question is whether the topographic features of different catchments will lead to distinct performance of the meshes generated by such LMR that is whether different refinement will be observed in relatively flat catchments compared to mountainous catchments with large topographic gradients.

Another key consideration in designing a suitable LMR strategy is to overcome inflexibility trade-offs (Ozdemir *et al.* 2013); these include coping with the need for a number of manually defined parameters to separately drive coarsening vs. refinement and with cross-parameter sensitivity issues and error control. Recent findings suggest that this type of inflexibility trade-off can be significantly reduced by taking a wavelet-based LMR (wLMR) strategy (Caviedes-Voullième & Kesserwani 2015). A wLMR strategy is founded on the theory of multiresolution analysis, which allows to systematically decide resolution levels by measuring a set of encoded details featuring the modeled data of interest relative to an error threshold chosen by the user. Starting from the finest available resolution, the multiresolution analysis gradually coarsens the resolution level until the evaluation of the encoded details indicates otherwise. Commonly, wavelet-based AMR (and therefore wLMR) are more suited for structured meshes (Cohen 2002). This limitation enables to recursively apply the same transforms to encode and decode the modeled data of interest across mesh resolution scales, hence making wLMR popular for quadtree-type meshes (for example, Gerhard *et al.* 2015). However, the use of wLMR to generate unstructured triangular meshes is yet to be explored, to find out how far it can favor multiresolution mesh generation for large-scale hydrological modeling applications.

In this context, triangular meshes still have certain desirable properties such as better capture complex domain and sub-catchment boundaries (Ivanov *et al.* 2004; Heinzer *et al.* 2012) and adaptability to match stream networks (Ivanov *et al.* 2004). LMR, to generate multiresolution unstructured triangular meshes (for example, Caviedes-Voullième *et al.* 2012; Ferraro *et al.* 2020), are relatively less common than structured quadtree-type meshes (for example, Heinzer *et al.* 2012; Hou *et al.* 2018) and could be more algorithmically complex to generate (Heinzer *et al.* 2012).

The first objective of this work is to propose and verify a strategy for applying wLMR to quad-based (that is structured grid-based) DEM data in order to generate multiresolution unstructured triangular meshes (hereinafter referred to as wLMR meshes) for catchment-scale hydrological simulations. The second objective is to carry out a comparative study of different input types for the criteria for wLMR. The performance of meshes generated by wLMR using topographic slope and topographic curvature as input is compared for rainfall–runoff events in a flat catchment (Thies catchment, Senegal) and in a mountainous catchment (Lower Triangle catchment of the East River, Colorado, USA).

## WAVELET-BASED LMR

The starting point is a rectangular domain Ω of length *L*, and width *W*, where the topography is given by uniform resolution DEM data. Initially, Ω is discretized by a fine resolution mesh *M* at the highest resolution available for the DEM data. From this fine resolution mesh, wLMR is applied to construct a multiresolution quadtree-type mesh *M _{s}* that has a reduced number of elements compared to

*M*. In this work,

*M*is constructed such that the perturbation error representing the deviation relative to

_{s}*M*(Gerhard

*et al.*2015) is forced to remain below a user-specified error threshold

*ɛ*. Performing multiresolution analysis to predict

*M*has been applied in a 1D × 1D manner using the Haar wavelet by following the procedure reported in Kesserwani

_{s}*et al.*(2019).

## PREDICTION OF *M*_{s}

_{s}

Starting from *M*, the maximum number of refinement levels *L*_{max} is assigned in order to perform multiresolution analysis to predict the refinement levels of *M _{s}*. A suitable choice of

*ɛ*is accepted to be proportional to the finest feature that needs capturing in the modeling (Gerhard

*et al.*2015; Kesserwani

*et al.*2019), which is here the finest resolution on

*M*. The coarsest length of mesh resolution on

*M*relative to Ω is denoted by Δ

_{s}*x*, which is chosen such that Δ

*x*

*=*

*W*if

*W*

*<*

*L*, or otherwise, it is set to Δ

*x*=

*L.*A cell size on the mesh

*M*is related to the coarsest mesh size Δ

*x*by

*L*

_{max}dyadic subdivision. On the multiresolution mesh

*M*, cell size can be at refinement level (

_{s}*n*), , with level (0) denoting the coarsest mesh possible. In this work,

*L*

_{max}is chosen such that the cell size at the refinement level (

*L*

_{max}) equals the available DEM data resolution.

*M*is represented as an

*r*×

*c*matrix

*A*, the discrete Haar wavelet transform can be recursively applied to produce the encoded differences between the DEM data across the two subsequent levels (

*n*

*+*1) and (

*n*) (e.g., Ruch & van Fleet 2009):

*W*in Equation (1) denotes the dimension of the matrix. That is to say that given the dimension

*d*,

*W*is a

_{d}*d*×

*d*matrix, consisting of the low-pass filter matrix

*H*and the high-pass filter matrix

_{d}*G*aswhere

_{d}*H*and

_{d}*G*are both

_{d}*d*/2 ×

*d*matrices. The transform

*B*has the structurewhere all , , , and are

*r*/2 ×

*c*/2 matrices. contains the

*scale coefficients*, whereas , , and contain the

*detail coefficients*(or details) in the vertical, horizontal, and diagonal directions, respectively. These detail coefficients are denoted by

*δ*

_{ij}^{(n)}and are analyzed to decide

*M*considering the details across all the resolution levels (

_{s}*n*),

*n*

*<*

*L*

_{max}. Here,

*i*and

*j*are the local row and column indices specifying a certain entry in all the matrices , , , and . Analysis of detail coefficients is performed in terms of their normalized magnitude (Gerhard

*et al.*2015):and in a level-wise manner starting from the resolution level (0) to (

*L*

_{max}− 1). Meanwhile, a detail coefficient is flagged as significant if its normalized detail satisfies:

Detailed explanation on how the resolution selection works in the context of a wavelet-based approach can be found in Kesserwani *et al.* (2019, Sec. 2.3.1).

## MULTIRESOLUTION TRIANGULAR MESH GENERATION

The quadtree-type multiresolution mesh *M _{s}* produced by the multiresolution analysis is used to constrain and trigger refinement to generate a multiresolution triangular mesh. The workflow is illustrated in Figure 1. Starting from a uniform triangular mesh at the coarsest refinement level (that is (

*n*) = (0)), each element of

*M*is mapped to a triangular cell. Each triangular cell containing an element of

_{s}*M*with a finer predicted resolution level than itself is flagged for refinement (indicated in Figure 1 by a dotted line). The flagged cells are refined once, and the process is repeated until no further triangular cell is flagged for refinement. This results in an unstructured triangular mesh version of

_{s}*M*.

_{s}## SURFACE FLOW NUMERICAL SOLVER

*M*, the surface flow is modeled by solving the two-dimensional zero-inertia equation on an unstructured triangular mesh. The equation writes where

_{s}*h*is the (ponding) water depth,

*t*denotes the time axis, ∇ is the Laplacian,

*U*is the surface flow velocity,

_{ZI}*s*is the rainfall source term, Γ

_{r}*is the water exchange between surface and subsurface systems,*

_{e}*n*is Manning's roughness coefficient,

*z*is the bed elevation, ||∇

*z*|| is the Euclidean

*L*-norm of ∇

*z*, and

*η*

*=*

*h*+

*z*is the free surface elevation of the water. In this work, the Advanced Terrestrial Simulator (ATS) (Coon

*et al.*2016) is used to numerically solve Equation (8) via an implicit finite-volume scheme.

## DIAGNOSTIC INVESTIGATION

This investigation aims to diagnostically explore the applicability of the proposed wLMR strategy for mesh refinement while comparing two types of DEM-based input data for the multiresolution analysis: elevation data, which will be interpreted as a slope-based refinement criterion for wLMR; and, bed slope data, which will be interpreted as a curvature-based refinement criterion for wLMR. Hereinafter, error thresholds for slope-based and curvature-based criteria are denoted with *ɛ ^{s}* and

*ɛ*, respectively.

^{c}The parameters used to set up and run these simulation runs are summarized in Table 1. Hydrographs at the outlet of the domain obtained on multiresolution meshes resulting from wLMR (wLMR meshes) are compared with reference hydrographs obtained from an alternative simulation on uniform meshes at the resolution Δ*x*^{(Lmax)}, that is to say the resolution of the highest refinement level. These results are compared using the root mean square error (RMSE) and the Nash–Sutcliffe model efficiency (NSE) coefficient (Nash & Sutcliffe 1970).

Symbol . | Parameter . | Value . |
---|---|---|

L | Length | 22 m |

Z | Bed elevation | See Figure 2 (left) |

C | Chézy coefficient | 1.767 ms^{−1/2} |

i | Rainfall intensity | See Figure 2 (top, right) |

ɛ, ^{s}ɛ ^{c} | Error thresholds | 0.01 m/m, 0.01 m^{−1} |

T | Time horizon | 4,800 s |

Γ(0) | Boundary condition at x= 0 | Closed |

Γ(L) | Boundary condition at x=L | Critical depth |

Symbol . | Parameter . | Value . |
---|---|---|

L | Length | 22 m |

Z | Bed elevation | See Figure 2 (left) |

C | Chézy coefficient | 1.767 ms^{−1/2} |

i | Rainfall intensity | See Figure 2 (top, right) |

ɛ, ^{s}ɛ ^{c} | Error thresholds | 0.01 m/m, 0.01 m^{−1} |

T | Time horizon | 4,800 s |

Γ(0) | Boundary condition at x= 0 | Closed |

Γ(L) | Boundary condition at x=L | Critical depth |

The simulation cases are chosen to be based on the analytical test reported in Govindaraju *et al.* (1988), which is characterized by a constant topographic slope and the hyetograph shown in Figure 2 (black line at the top right). In all the simulations, *ɛ ^{s}* and

*ɛ*are both set to be in the order of magnitude of the cell size corresponding to the finest resolution (around

^{c}*ɛ*= 0.01 m/m;

^{s}*ɛ*= 0.01 m

^{c}^{−1}).

Predicted refinement levels for all cases are plotted in Figure 2 (dashed and dotted lines at the top right). For the constant topographic slope, the slope-based criterion predicts uniform refinement at level (6), just one level below *L*_{max} = (7). Similar behavior is reported in Kesserwani *et al.* (2019) for a dam break over sloped bed. In contrast, the curvature-based refinement criterion predicts no refinement at all, which implies that in the absence of curvature, the curvature-based refinement criterion is not fit for purpose. In comparison, for constant topographic curvature, the slope-based refinement criterion suggests a coarser refinement level (3) at the beginning of the domain where the slope is small, and increases the refinement level up to (6) at the outlet of the domain where the slope is steep, while the curvature-based refinement criterion predicts a constant refinement level (6) (similar to the slope-based refinement criterion applied to the constant topographic slope). For the variable topographic curvature, both topographic slope and curvature are steep at the beginning of the domain and small at the outlet. Thus, both slope-based and curvature-based refinement criteria predict a maximum level of refinement at the beginning of the domain.

Computed hydrographs at the outlet of the domain are plotted in Figure 2 (right) and are summarized in Table 2. In all cases, the model results do not significantly deviate from each other and show good agreement with the reference solution. In the case with constant slope, the curvature-based criterion yields a poor fit to the reference hydrograph due to numerical diffusion caused by the coarse resolution.

Case . | RMSE (m^{2}s^{−1})
. | NSE . | N (cells)
. |
---|---|---|---|

CS (s) | 9.0 × 10^{−3} | 0.99 | 512 |

CS (c) | 9.4 × 10^{−2} | 0.99 | 8 |

CC (s) | 6.4 × 10^{−4} | 0.99 | 416 |

CC (c) | 5.4 × 10^{−4} | 0.99 | 512 |

VC (s) | 3.2 × 10^{−3} | 0.99 | 428 |

VC (c) | 6.5 × 10^{−3} | 0.99 | 366 |

Case . | RMSE (m^{2}s^{−1})
. | NSE . | N (cells)
. |
---|---|---|---|

CS (s) | 9.0 × 10^{−3} | 0.99 | 512 |

CS (c) | 9.4 × 10^{−2} | 0.99 | 8 |

CC (s) | 6.4 × 10^{−4} | 0.99 | 416 |

CC (c) | 5.4 × 10^{−4} | 0.99 | 512 |

VC (s) | 3.2 × 10^{−3} | 0.99 | 428 |

VC (c) | 6.5 × 10^{−3} | 0.99 | 366 |

The letter inside the parentheses denotes the steering parameter (i.e., s: slope, c: curvature).

These test cases confirm that the wLMR strategy behaves as expected for different input data. Using the bed elevation as input for the wLMR strategy results in refinement based on the slope, while using the slope tensor as input results in refinement based on the curvature. Further, it is observed that neither criterion clearly outperforms the other. The performance of the refinement criteria depends on the characteristics of the topography. In the absence of curvature, using the curvature-based refinement criterion is not advisable. This suggests that the shape of the topography must be considered when choosing the type of refinement criterion.

## CASE STUDIES AT THE CATCHMENT SCALE

These case studies explore and contrast the applicability of the proposed wLMR strategy to different real-world topography. In particular, the performance of the wLMR strategy with its two refinement criteria is compared for two catchments with distinct features: the Thies catchment in Senegal, which is a laboratory-scale catchment characterized by smooth topography; and the Lower Triangle catchment in Colorado, USA, which is characterized by rugged mountainous topography. The wLMR is applied with a range of values for *ɛ ^{s}* and

*ɛ*, chosen such that their resulting meshes have approximately the same number of cells, denoted by

^{c}*N*. Reference solutions are obtained on uniform triangular meshes at available data resolution. On each mesh, the predicted topography and velocity profiles are evaluated relative to the reference solutions using quantitative indices.

## THIES CATCHMENT, SENEGAL

This case considers steady-state flow in a laboratory-scale catchment in Thies, Senegal (Mügler *et al.* 2011). Figure 3(a) shows the topography of the catchment, based on DEM data with 0.1 m resolution. The catchment is characterized by smooth topography; steep slopes are only present at the outlet of the domain. Measurements of steady-state flow velocities are available at 62 gauges inside the catchment. The black dots in Figure 3(b) (top) show the locations of these gauges.

Table 3 summarizes the parameters used to set up and run the simulation, taken from Mügler *et al.* (2011). The reference solution is generated by a simulation run on an unstructured mesh with a uniform resolution of 0.1 m (*N* = 7,888 cells). Figure 3(b) (top) shows the magnitude of these reference velocities. The flow concentrates in the middle of the domain, where a flat channel is fed by several smaller streams from the side. Figure 3(b) (bottom) plots the correlation between reference velocities and measurements. The reference model predicts the velocities reasonably well; the RMSE is about 0.038 m s^{−1}. These results are consistent with the results in Simons *et al.* (2014) and confirm the validity of the reference solution.

Symbol . | Parameter . | Value . |
---|---|---|

L | Length | 10 m |

W | Width | 4 m |

z | Bed elevation | see Figure 3(a) |

n | Manning coefficient | 0.02 ms^{−1/2} |

i | Rainfall intensity | 51.5 × 10^{−3} m h^{−1} |

ɛ ^{s}, ɛ^{c} | Error threshold | See the text |

T | Time horizon | Until steady state |

Γ(0) | Boundary condition at x= 0 | Closed |

Γ(L) | Boundary condition at x=L | Critical depth |

Symbol . | Parameter . | Value . |
---|---|---|

L | Length | 10 m |

W | Width | 4 m |

z | Bed elevation | see Figure 3(a) |

n | Manning coefficient | 0.02 ms^{−1/2} |

i | Rainfall intensity | 51.5 × 10^{−3} m h^{−1} |

ɛ ^{s}, ɛ^{c} | Error threshold | See the text |

T | Time horizon | Until steady state |

Γ(0) | Boundary condition at x= 0 | Closed |

Γ(L) | Boundary condition at x=L | Critical depth |

We then generate wLMR meshes using the different thresholds listed in Table 4. Figure 4 shows the predicted refinement levels for both the slope- and curvature-based refinement criteria needed to generate the same amount of cells *N*. Figure 4(a) shows the refinement levels predicted by the slope-based refinement criterion for *N* = 300, 450, 600, and 2,300 cells; whereas Figure 4(b) shows those predicted by the curvature-based refinement criterion for these same *N*. Both refinement criteria predict different distributions of refinement levels: in Figure 4(a), we see that the slope-based refinement criterion refines in the steep regions. The channel in the middle of the domain is coarsened, because the topography in this region is flat and smooth. In comparison, Figure 4(b) shows that the curvature-based refinement criterion refines around the channel and its contributing streams. The slope-based refinement criterion coarsens more aggressively than the curvature-based criterion. However, the refinement of the curvature-based criterion is more sensible, because it refines around the channel where the flow concentrates.

ɛ (m/m)
. ^{s} | N (cells)
. ^{s} | RMSE elevation (m)
. ^{s} | RMSE slope (m/m)
. ^{s} | RMSE velocity (ms^{s}^{−1})
. | PRF . |
---|---|---|---|---|---|

(a) | |||||

1 × 10^{−6} | 2,437 | 0.0013 | 2.5 × 10^{−4} | 0.017 | 2.5 |

5 × 10^{−6} | 584 | 0.0022 | 3.7 × 10^{−4} | 0.018 | 2.8 |

1 × 10^{−5} | 450 | 0.0026 | 4.4 × 10^{−4} | 0.020 | 3.0 |

5 × 10^{−5} | 300 | 0.0045 | 6.2 × 10^{−4} | 0.021 | 3.6 |

ɛ (m^{c}^{−1})
. | N (cells)
. ^{c} | RMSE elevation (m)
. ^{c} | RMSE slope (m/m)
. ^{c} | RMSE velocity (ms^{c}^{−1})
. | PRF . |

(b) | |||||

1 × 10^{−6} | 2,270 | 0.0011 | 2.2 × 10^{−4} | 0.012 | 2.6 |

5 × 10^{−6} | 606 | 0.0025 | 3.9 × 10^{−4} | 0.015 | 2.8 |

1 × 10^{−5} | 445 | 0.0028 | 4.1 × 10^{−4} | 0.016 | 3.0 |

5 × 10^{−5} | 320 | 0.0033 | 5.3 × 10^{−4} | 0.018 | 3.6 |

ɛ (m/m)
. ^{s} | N (cells)
. ^{s} | RMSE elevation (m)
. ^{s} | RMSE slope (m/m)
. ^{s} | RMSE velocity (ms^{s}^{−1})
. | PRF . |
---|---|---|---|---|---|

(a) | |||||

1 × 10^{−6} | 2,437 | 0.0013 | 2.5 × 10^{−4} | 0.017 | 2.5 |

5 × 10^{−6} | 584 | 0.0022 | 3.7 × 10^{−4} | 0.018 | 2.8 |

1 × 10^{−5} | 450 | 0.0026 | 4.4 × 10^{−4} | 0.020 | 3.0 |

5 × 10^{−5} | 300 | 0.0045 | 6.2 × 10^{−4} | 0.021 | 3.6 |

ɛ (m^{c}^{−1})
. | N (cells)
. ^{c} | RMSE elevation (m)
. ^{c} | RMSE slope (m/m)
. ^{c} | RMSE velocity (ms^{c}^{−1})
. | PRF . |

(b) | |||||

1 × 10^{−6} | 2,270 | 0.0011 | 2.2 × 10^{−4} | 0.012 | 2.6 |

5 × 10^{−6} | 606 | 0.0025 | 3.9 × 10^{−4} | 0.015 | 2.8 |

1 × 10^{−5} | 445 | 0.0028 | 4.1 × 10^{−4} | 0.016 | 3.0 |

5 × 10^{−5} | 320 | 0.0033 | 5.3 × 10^{−4} | 0.018 | 3.6 |

We also evaluate the predicted topography to determine possible a priori indicators for the accuracy of the terrain connectivity prediction on wLMR meshes. First, we evaluate the significance of bed elevations. Figure 5 shows the correlation between bed elevations predicted on the wLMR meshes and the uniform reference mesh. It plots the bed elevation of each cell in the wLMR mesh against its corresponding reference bed elevation. For a perfect correlation, all points would be located on the black line. Points with a long distance to this black line are poorly approximated. Elevations obtained on meshes generated by the slope-based criterion are plotted with circles; elevations obtained on meshes generated by the curvature-based criterion are plotted with triangles. Although both criteria refine in different regions (recall Figure 4), they both accurately reproduce the reference topography for all *N*. Thus, bed elevation might not be a significant indicator for the accuracy of the terrain connectivity prediction on wLMR meshes.

Another possible a priori indicator for the accuracy of the terrain connectivity prediction on wLMR meshes is the bed slope. Figure 6 shows the correlation between the bed slopes predicted on wLMR and reference meshes. The deviation between these bed slopes is significant, but improves with increasing *N*; see the RMSE in Table 4. Further, in contrast to bed elevation, both refinement criteria show different behaviors for slope. In the bottom part of Figure 6, we see that for *N* > 600 cells, the bed slopes predicted on the mesh generated by the curvature-based refinement criterion are steeper than the bed slopes predicted on the mesh generated by the slope-based criterion. Thus, the bed slope is likely a more reliable indicator to measure the capability of the wLMR to capture channeling patterns.

We now run the simulations on wLMR meshes and compare the steady-state velocities on the wLMR mesh relative to the reference solution at the 62 gauges. Figure 7 shows the correlation between these velocities. The velocities agree satisfactorily with each other. Same as bed elevation and slope, the accuracy of the predicted velocities increases with increasing *N*; see the RMSE in Table 4, which also shows that the curvature-based refinement criterion performs better than the slope-based wLMR.

Evaluating the predicted bed slopes indicates that the curvature-based refinement criterion captures the topography characteristics more accurately; see Figure 6, which results in a more accurate prediction of flow velocities inside the domain. This suggests that the accuracy of the bed slope prediction can be used as an a priori indicator for the accuracy terrain connectivity prediction on wLMR meshes. Specifically, evaluating the predicted bed slopes may guide the choice of input type – slope or curvature – for the refinement criteria of the wLMR.

In this case, for both refinement criteria, the approximation accuracy for all evaluated values improves with increasing *N*. This is because the available elevation data have a sufficient high resolution, which enables the refinement criteria to detect small-scale topography features and successively refine the mesh in these areas.

Table 4 also lists the performance reduction factor (PRF), which is the ratio between the CPU times of the high-resolution reference simulation and the simulations run on the wLMR meshes. A large value of PRF indicates that the simulations run on the wLMR meshes reduce the CPU time significantly. The observed reduction in CPU time is not overwhelming.

## LOWER TRIANGLE CATCHMENT OF THE EAST RIVER, COLORADO, USA

In this case study, transient flow is considered in the Lower Triangle catchment of the East River, CO, USA (Carroll *et al.* 2018; Hubbard *et al.* 2018) for which a DEM is available at 10 m resolution. This is a mountainous catchment with rugged topography and sudden and large local variations in slopes, in contrast with the flat Thies catchment studied in the previous section that also involved a DEM at a much finer resolution. Figure 8(a) shows the topography of the catchment that has a surface area of 14.82 km^{2}, with a bed elevation ranging from 2,759 m at its south-western part to 3,787 m at its northern part. At lower elevations – darker shades of gray in Figure 8(a) – a flat river channel is located in the middle of the catchment. The refinement criterion needs to capture this river channel for an accurate prediction of flow inside the domain.

Simulations were run on a uniform unstructured mesh and several wLMR meshes were obtained by slope- and curvature-based refinement criteria, using the setup parameters listed in Table 5. The rainfall data used in these simulations correspond to a high-intensity rainfall event extracted from a 1-year rainfall time series with 30 min resolution; see Figure 8(b). The setup parameters and rainfall follow (Dwivedi *et al.* 2018). Results obtained on the uniform unstructured mesh (cell size of 10 m; *N* = 291,614 cells) are considered as the reference results. The deviation between the predictions obtained on the wLMR meshes and the reference results is studied as a means to compare slope- and curvature-based refinement criteria. Due to the lack of available discharge measurement data, the model could not be calibrated. However, an accurate calibration is not strictly necessary within the scope of comparing predictions obtained on wLMR meshes.

Symbol . | Parameter . | Value . |
---|---|---|

z | Bed elevation | See Figure 8(a) |

n | Manning coefficient | 0.03 ms^{−1/2} |

i | Rainfall intensity | See Figure 8(b) |

ɛ | Error threshold | Various, see text |

T | Time horizon | 1.5 days |

Γ(z < 2,700 m) | Boundary condition at the outlet | Critical depth |

Γ(else) | Boundary condition on all other sides | Closed |

Symbol . | Parameter . | Value . |
---|---|---|

z | Bed elevation | See Figure 8(a) |

n | Manning coefficient | 0.03 ms^{−1/2} |

i | Rainfall intensity | See Figure 8(b) |

ɛ | Error threshold | Various, see text |

T | Time horizon | 1.5 days |

Γ(z < 2,700 m) | Boundary condition at the outlet | Critical depth |

Γ(else) | Boundary condition on all other sides | Closed |

For the meshes generated by the two refinement criteria, the thresholds listed in Table 6 were chosen such that the slope-based and curvature-based refinement criteria result in meshes with the same number of cells *N*, also given in Table 6. For these *N*, Figure 9(a) shows the refinement levels predicted by the slope-based refinement criterion; whereas Figure 9(b) shows those predicted by the curvature-based refinement criterion. In terms of coarsening, the curvature-based approach coarsens more aggressively, but neither of the approaches could capture the channel in the middle of the domain at the finest resolution. Analysis of refinement levels does not reveal whether the slope-based or the curvature-based refinement criteria give a more sensible refinement for this test.

ɛ (m/m)
. ^{s} | N (cells)
. ^{s} | RMSE elevation (m)
. ^{s} | RMSE slope (m/m)
. ^{s} | RMSE discharge (m^{s}^{3}s^{−1})
. | NSE
. ^{s} | PRF . |
---|---|---|---|---|---|---|

(a) | ||||||

7.5 × 10^{−4} | 190,024 | 2.50 | 0.1 | 124.5 | 0.93 | 1.9 |

1.0 × 10^{−3} | 161,568 | 1.72 | 0.09 | 198.7 | 0.81 | 2.1 |

1.8 × 10^{−3} | 95,191 | 1.42 | 0.08 | 210.6 | 0.80 | 2.6 |

5.0 × 10^{−3} | 42,449 | 1.35 | 0.08 | 264.7 | 0.67 | 3.2 |

ɛ (m^{c}^{−1})
. | N (cells)
. ^{c} | RMSE elevation (m)
. ^{c} | RMSE slope (m/m)
. ^{c} | RMSE discharge (m^{c}^{3}s^{−1})
. | NSE
. ^{c} | PRF . |

(b) | ||||||

0.1 | 183,890 | 2.63 | 0.1 | 96.4 | 0.96 | 2.0 |

0.125 | 164,751 | 1.82 | 0.08 | 110.0 | 0.95 | 2.2 |

0.25 | 98,189 | 1.47 | 0.08 | 122.2 | 0.94 | 2.6 |

0.5 | 47,653 | 1.40 | 0.08 | 257.5 | 0.72 | 3.1 |

ɛ (m/m)
. ^{s} | N (cells)
. ^{s} | RMSE elevation (m)
. ^{s} | RMSE slope (m/m)
. ^{s} | RMSE discharge (m^{s}^{3}s^{−1})
. | NSE
. ^{s} | PRF . |
---|---|---|---|---|---|---|

(a) | ||||||

7.5 × 10^{−4} | 190,024 | 2.50 | 0.1 | 124.5 | 0.93 | 1.9 |

1.0 × 10^{−3} | 161,568 | 1.72 | 0.09 | 198.7 | 0.81 | 2.1 |

1.8 × 10^{−3} | 95,191 | 1.42 | 0.08 | 210.6 | 0.80 | 2.6 |

5.0 × 10^{−3} | 42,449 | 1.35 | 0.08 | 264.7 | 0.67 | 3.2 |

ɛ (m^{c}^{−1})
. | N (cells)
. ^{c} | RMSE elevation (m)
. ^{c} | RMSE slope (m/m)
. ^{c} | RMSE discharge (m^{c}^{3}s^{−1})
. | NSE
. ^{c} | PRF . |

(b) | ||||||

0.1 | 183,890 | 2.63 | 0.1 | 96.4 | 0.96 | 2.0 |

0.125 | 164,751 | 1.82 | 0.08 | 110.0 | 0.95 | 2.2 |

0.25 | 98,189 | 1.47 | 0.08 | 122.2 | 0.94 | 2.6 |

0.5 | 47,653 | 1.40 | 0.08 | 257.5 | 0.72 | 3.1 |

The RMSE for discharge and the NSE are calculated based on the hydrograph at the outlet of the domain.

Similar to the previous case, comparing the deviations between the bed elevation quantity on both wLMR meshes relative to that on the reference mesh indicates that both refinement criteria accurately capture the reference bed elevation for all *N*. As the differences between the bed elevations predicted on the wLMR meshes were visually indistinguishable, they were not illustrated and studied. This can be reinforced by comparing the RMSE in Table 6 for the bed elevations on both wLMR meshes. Again, comparing bed elevations is found not to be a significant a priori indicator for measuring the accuracy of the terrain connectivity prediction on a wLMR mesh. Therefore, the correlation between model predictions on wLMR and reference meshes is further analyzed in terms of bed slopes – shown in Figure 10 – to assess whether the bed slope is a better a priori indicator. As shown in Figure 10, both refinement criteria now predict very similar bed slopes that nonetheless significantly deviate from the reference bed slopes. It can also be noted that increasing *N* gives a marginal improvement in the agreement between results obtained on wLMR and reference meshes. Bed slopes beyond *N* = 96,000 cells do not improve significantly. This is also seen in the RMSE for bed slopes in Table 6.

Figure 11 compares the transient hydrographs at the outlet of the catchment computed on the slope- and curvature-based wLMR meshes and the uniform reference mesh. The hydrographs predicted on the wLMR meshes are generally in a good agreement with the reference hydrograph. The hydrograph predicted on the mesh generated by the curvature-based refinement criterion more accurately predicts the first peak observed in the reference hydrograph at around 5 h, whereas the hydrograph on the mesh generated by the slope-based refinement criterion falls short. Both refinement criteria correctly capture the second peak in their predicted hydrographs at around 10 h, which is caused by the peak in the rainfall intensity. Similar to the bed slope prediction, the hydrograph predictions on both wLMR meshes do not significantly improve beyond *N* = 96,000 cells. The only significant improvement is observed in the hydrograph predicted on the wLMR mesh generated by the slope-based criterion, which improves its prediction of the first peak for *N* = 185,000 cells. The reduced improvement in accuracy beyond *N* = 96,000 cells is also quantitatively notable in the RMSE and NSE of the hydrographs obtained on wLMR meshes with regard to the reference model hydrograph in Table 6.

Figure 12 shows the correlation between velocities on the wLMR meshes and reference velocities on the uniform mesh at different output times. Despite the good agreement between the hydrographs on the two wLMR meshes, seen in Figure 11, their predicted velocities significantly deviate from reference velocities at all output times and for all *N*. At *t* = 5 h, the first peak of the hydrograph is reached. Overall, the velocities predicted on meshes generated by both the slope- and curvature-based refinement criteria show similar agreement with the reference velocities, see Figure 12(a). At *t* = 10 h, the second peak of the hydrograph is reached. As seen in Figure 12(b), at this time, the velocity predictions on meshes generated by the curvature-based refinement criterion show a better agreement with reference velocities on the uniform mesh than the velocity predictions on meshes generated by the slope-based refinement criterion. After *t* = 20 h, the rainfall has ceased and the catchment is draining the remaining water. The velocity predictions on meshes generated by the curvature-based refinement criterion grossly overestimate the reference velocities, while the velocity prediction on meshes generated by the slope-based refinement criterion shows a better agreement; see Figure 12(c). Table 6 summarizes the RMSE for velocities predicted on wLMR meshes with regard to the reference velocities for all *t* and *N*. We observe that velocities based on both refinement criteria have similar RMSE, with exception of the severe outliers of the curvature-based refinement criterion that results in gross overestimations of the reference velocities; see, for example, Figure 12(c). Comparison of the RMSE in Table 6 further indicates that increasing *N* beyond 96,000 cells does not significantly enhance the accuracy of the velocity predictions. However, depending on the value of *N*, velocity predictions on the curvature-based wLMR meshes seem to be prone to severe overestimations, such that going beyond 96,000 cells may worsen their agreement with reference velocities. In contrast, the velocity predictions on the meshes generated by the slope-based refinement criterion do not suffer from such overestimations.

The hydrograph prediction on both wLMR meshes shows good agreement with the reference hydrograph. But the good agreement between the hydrographs does not translate to the prediction of velocity – recall Figures 11 and 12. The hydrograph is a domain-integrated signature and is, in general, not sensitive to local spatial variability inside the domain, especially for high-flow regimes (Khosh Bin Ghomash *et al.* 2019). This means that even if local flow states inside the domain are predicted poorly, the hydrograph might still be predicted well.

On the other hand, the accurate prediction of local flow states, such as the flow velocities, requires capturing the flow connectivity of the domain, especially for low-flow regimes (Caviedes-Voullième *et al.* 2020). We see, in Figure 12, that local flow states are more robust to mesh resolution when the domain is fully submerged and terrain connectivity becomes less significant, which is the case during the second peak in Figure 12(b). In contrast, terrain connectivity becomes significant when the catchment is drained and thus, velocities in Figure 12(c) are sensitive to mesh resolution. This is most likely due to the representation of terrain connectivity on the wLMR meshes (Caviedes-Voullième *et al.* 2020). The (lack of) accuracy in the predicted terrain connectivity also may explain the improvement in the predicted hydrographs for *N* = 185,000 cells in Figure 11. Increasing *N* leads to additional refinement in regions near to the outlet, which enhances the representation of the terrain connectivity in these regions. Thus, the first peak of the reference hydrograph is captured more accurately, which also leads to a better prediction of the second peak – especially on the slope-based wLMR meshes.

For increasing *N*, a rapidly diminishing improvement in the prediction accuracy on both wLMR meshes is observed for all predicted quantities. A reason for this might be the constraint posed by the available topography data resolution of 10 m. This resolution is too coarse to capture significant small-scale variations in flat regions. Hence, lowering the threshold does not lead to refinement in these areas. We saw similar behavior in the diagnostic investigation, where in the absence of curvature, the curvature-based refinement criterion predicted no refinement. This means that, depending on the resolution and characteristics of the available topography data, there is a certain threshold beyond which no significant improvement can be gained. The value of this threshold can be estimated by evaluating the accuracy of the bed slope prediction on the wLMR meshes; see Figure 10.

The PRF for all *N* is shown in Table 6. As observed in the previous case, the reduction in CPU time is again low; it ranges between 2 and 3.

## CONCLUSIONS

This work adapted a wavelet-based LMR (wLMR) strategy to generate multiresolution meshes with unstructured triangular elements from quad-based DEM for rainfall–runoff simulations. The wLMR strategy used as input either the bed elevation of the DEM or its bed slopes that led to two refinement criteria: a slope-based refinement criterion and a curvature-based refinement criterion, respectively.

The numerical model used in this work (ATS) simulated the overland flow dynamics by solving the two-dimensional zero-inertia equation – see Equation (8) – with an implicit finite-volume scheme. Simulations of rainfall–runoff were carried out for idealized geometries and two realistic catchments with contrasting flow and topography features – steady state vs. transient; flat and smooth vs. steep and rugged – were carried out. The model results obtained on meshes generated by the wLMR strategy for both criteria were systematically compared to results obtained on uniform resolution meshes at the available data resolution. The bed elevation, bed slope, and flow velocity predictions were compared in a cell-by-cell manner using RMSE as a quantitative metric. Plots of the correlation between predictions on different meshes were also provided for the discussion. For cases featuring transient flow, the hydrograph at the outlet of the domain predicted on wLMR meshes was compared to a reference hydrograph obtained on a high-resolution uniform mesh in terms of RMSE and NSE. While flow velocity and hydrograph predictions were used to evaluate the performance of the wLMR meshes in terms of runoff prediction capability, bed elevation and bed slope predictions were used to assess the suitability of these terrain characteristics as a priori indicators for the accuracy of the terrain connectivity prediction on wLMR meshes.

Simulation runs on idealized geometries showed that the slope- and curvature-based refinement criteria are equally effective and that the shape of the topography significantly affects the performance of the wLMR strategy. Where the topography has no curvature – a constant slope, the curvature-based refinement criterion is unable to detect critical areas and should not be used. In these cases, using a slope-based wLMR gives a more sensible refinement.

For the real case studies at laboratory- and real-scale, the bed elevation was well predicted on meshes generated by both refinement criteria for all values of *N*. In contrast, the bed slope was not predicted to the same extent, which might reflect the scale-dependency of the topographic slope. In these cases, bed slope analysis provided a better metric of the correlation between a wLMR mesh and a uniform mesh. While the hydrograph at the outlet of the domain was well predicted on all wLMR meshes, velocity predictions inside the domain deviated from the high-resolution reference velocities. This suggests that the terrain connectivity is not reproduced correctly by any of the wLMR meshes. In the flat and smooth catchment, the curvature-based refinement criterion was found to better capture the channeling patterns, while in the steep and rugged catchment, neither of the refinement criteria clearly captured these. Nevertheless, the curvature-based refinement criterion predicted the hydrograph better, while the slope-based refinement criterion predicted the velocities inside the domain more consistently. Overall, the performance of the refinement criteria was found to be dependent on the shape and resolution of the available topography data.

In the investigated cases, using wLMR on the same number of CPUs as the high-resolution reference simulation reduces computational cost roughly by a factor between 2 and 3. This rather small increase in computational efficiency is likely due to the simplified governing equations, communication overhead of the distributed parallelization, and the use of an implicit solver. A greater efficiency from using wLMR may be more apparent in more computationally expensive simulations, such as in domains with millions of cells or in problems considering coupled processes. For example, if overland flow is coupled to subsurface hydrology – as is increasingly needed in science applications (Hubbard *et al.* 2018) – the resulting problems are much larger. It is foreseeable that the advantage of mesh refinement in such problems may become more significant. For rainfall-induced surface runoff processes, as considered in this work, the terrain connectivity was found to be an important property to be preserved by the AMR. In other types of problems such as river flooding or coupled subsurface flow, different properties may need to be preserved by the AMR strategy. As such, it is likely that different conclusions than in this work will be reached. Thus, the applicability of the presented wLMR strategy remains to be explored in a range of problems beyond the one presented here. However, a key aspect of this strategy is that it is flexible and general enough that other data may be used as input to generate multiresolution meshes for different objectives.

## ACKNOWLEDGEMENTS

This work is funded as part of IDEAS Watersheds by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research (award no. DE-AC02-05CH11231). G. Kesserwani acknowledges the UK Engineering and Physical Sciences Research council (grant ID: EP/R007349/1). The simulations of the Lower Triangle catchment of the East River, CO, USA, are based upon work supported as part of the Sustainable Systems Scientific Focus Area funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research (award no. DE-AC02-05CH11231). This work used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy, Office of Science, user facility operated under contract no. DE-AC02-05CH11231.