Joint Effects of the DEM Resolution and the Computational Cell Size on the Routing Methods in Hydrological Modelling

Abstract

Natural disasters, including droughts and floods, have caused huge losses to mankind. Hydrological modelling is an indispensable tool for obtaining a better understanding of hydrological processes. The DEM-based routing methods, which are widely used in the distributed hydrological models, are sensitive to both the DEM resolution and the computational cell size. Too little work has been devoted to the joint effects of DEM resolution and computational cell size on the routing methods. This study aims to study the joint effects of those two factors on discharge simulation performance with two representative routing methods. The selected methods are the improved aggregated network-response function routing method (I-NRF) and the Liner-reservoir-routing method (LRR). Those two routing methods are combined with two runoff generation models to simulate the discharge. The discharge simulation performance is evaluated under the cross combination of four DEM resolutions (i.e., 90 m, 250 m, 500 m, and 1000 m) and fifty-six computational cell sizes (ranging from 5 arc-min to 60 arc-min). Eleven years of hydroclimatic data from the Jianxi basin (2000–2010) and the Shizhenjie basin (1983–1993) in China are used. The results show that the effects of the DEM resolution and the computational cell size are different on the I-NRF method and the LRR method. The computational cell size has nearly no influence on the performance of the I-NRF methods, while the DEM resolution does. On the contrary, the LRR discharge simulation performance decreases with oscillating values as the computational cell size increases, but is hardly affected by the DEM resolution. Furthermore, the joint effects of the DEM resolution and the computational cell size can be ignored for both routing methods. The results of this study will help to establish the appropriate DEM resolution, computational cell size, and routing method when researchers build hydrological models to predict future disasters.

1. Introduction

Natural disasters bring huge losses to society, ecology, and the economy. Drought and flood events are common natural disasters and have led to significant troubles for humans [1,2,3]. Drought events can cause extreme shortages of domestic water, irrigation water and food, such as the drought that occurred in Southwest China in 2009 [4] and the droughts in eastern Africa in 2011 and 2017 [5], while extreme flood events have also caused casualties and property losses [6,7,8,9].

The development of science and technology [10,11,12,13,14,15] provides the possibility for disaster prevention and control. Studying the natural laws of the basin [16,17] and predicting the hydrometeorological elements of the basin are important measures to prevent natural disasters. Hydrological models play an important role in hydrological simulation and hydrological prediction. The accuracy of the hydrological simulation is affected by many factors, such as the routing methods [18], the accuracy of the input data [19,20], the model parameters [21] and so on [22]. The routing method, which transfers the simulated runoff over the basin into the downstream river discharge in hydrological modelling, plays a vital role in drought and flood prediction. A robust routing method contributes to achieve better discharge simulation results. The fluvial hydraulics model is a satisfactory method to establish the hydraulic hazard maps [23,24]. It simulates detailed flood dynamics at high spatial and temporal resolutions [24]. However, a lot of measured data, including hydrological data and hydraulic data, are needed in hydraulics simulation. Besides, highly sophisticated hydraulic simulations are time consuming [20], especially for large-scale watersheds. In the absence of detailed measured data or for basins of large areas, the DEM-based routing methods, or intelligent algorithms [25], can be alternatives. Conceptual models require less computational effort but they do not take into account the flood inundation and wave diffusion. The growing number of easily accessible gridded meteorological datasets and global Digital Elevation Model (DEM) datasets have facilitated the rapid development and widespread use of DEM-based routing methods. DEM-based routing methods are widely used in flood routing [23,26,27,28]. Samela et al. [23] proposed a DEM-based approach for extracting the flood-prone areas in the ungauged basin. It behaves well in a subcatchment of the Bulbula River in Africa. Li et al. [29] applied two DEM-based routing methods to the distributed Hydrologiska Byråns Vattenbalansavdelning (HBV) model. The discharge simulation performances with the DEM-based routing methods were better than that of the elevation band-based HBV model.

There are two types of DEM-based routing methods: Cell-to-Cell and Source-to-Sink [18,30]. Both are widely used in distributed hydrological models. The Cell-to-Cell method, which is based on both a channel storage equation and a water balance equation, determines the outflow of each cell, using the known inflow and storage at each time step. The outflow of one cell is the inflow of the adjacent downstream cell. These calculations are performed iteratively at each cell until the outflow of the basin is obtained. The main advantage of Cell-to-Cell routing methods is that the discharge is determined at every cell and the runoff is fully distributed over the basin. Cell-to-Cell routing methods are often fairly simple to implement and understand, although the computational demand increases rapidly with the increase in the number of cells [31,32,33,34]. The Source-to-Sink routing method determines how the runoff in each cell independently contributes to the total discharge at the outlet of the basin. Most Source-to-Sink routing methods are more computationally efficient than the Cell-to-Cell methods and are more widely applicable in various disciplines [30,35,36,37].

Generally, runoff is routed through flow nets identified from DEM data, in both the Cell-to-Cell and Source-to-Sink routing methods. Due to computational limitations in large-scale (i.e., large watershed, national or global) hydrological models, routing is often performed with computational cell resolutions that are coarser than those found in the DEMs.

In practice, the computational cell size may vary widely within a given model, depending on the spatial domain of the model works; while relatively coarser computational cell sizes may be used for simulating continental-scale basins, finer computational cell sizes are expected to improve the accuracy of the calibrated parameters and the model results. In regional and global hydrological modelling, the computational cell size usually ranges from 0.1° $\times$ 0.1° (about 10 km $\times 10 \mathrm{km}$) to 5$°\times$ 5$°$ (about 500 km $\times$ 500 km). For example, Harrigan et al. [38] simulated discharge at 1801 discharge stations, using the LISFLOOD routing model [39] and the ERA5 global runoff reanalysis dataset, which has a resolution of 0.1° × 0.1° [40]. Xie et al. [41] conducted a regional parameter estimation in China, using the Variable Infiltration Capacity (VIC) model at a resolution of 50 km $\times$ 50 km. Niu et al. [42] simulated the discharge based on the flow net at a resolution of 1$°\times$ 1$°$ in 16 catchments in the Pearl River Basin in China. Just as there are various computational cell size options, the typical DEM resolutions in hydrological modelling range from 15 m $\times$ 15 m to 1 km $\times$ 1 km. While finer DEM data resolutions are desirable, in practice, the preferred DEM resolution must be balanced against the available computational resources, especially in large-scale hydrological simulation. If the computational costs are too high, or if there is a lack of high-resolution DEM data for a given study area, it may make more sense to use a coarser DEM dataset.

The DEM dataset resolution and the computational cell size affect the routing method performance at the same time, and the user must carefully select the appropriate computational cell size and DEM resolution for their specific situation. While numerous studies have examined how the two conditions affect the modelling results separately [43,44], there is very little information about how they jointly affect the routing method performance. Using two different computational cell sizes of 25 km and 350 km, Arora et al. [45] explored how the computational cell size influenced the routing method results (one is for coupled general circulation models at the global scale and the other is for hydrological applications at basin scale). Overall, coarser computational cells resulted in poorer performances. Furthermore, Gong et al. [18] found that, in a water and snow balance modelling system (WASMOD) model of the Dongjiang Basin, the aggregated network-response function (NRF) performed similarly over a variety of computational cell sizes; conversely, the performance of the linear-reservoir-routing method (LRR) decreased with a zigzag pattern as the computational cell size increased. The DEM resolution influenced numerous hydrological parameters, such as the slope, the routing length, the flow net, and, consequently, the simulated discharge [46]. Using DEM resolutions of 100 m, 200 m, and 300 m, Du et al. [44] investigated the effects of the DEM resolution on the simulated storm runoff in a time variant spatially distributed travel time method, which showed that the performance was poor when the DEM resolution was larger than 200 m. Gong et al. [35] demonstrated that the simulated discharge, using the NRF routing method, was relatively insensitive to the DEM resolution differences between the Hydroshed (about 90 m) and Hydro1k (1000 m) datasets. Jeon et al. [47] studied the discharge results and processing time in the Heukcheon Basin for various DEM resolutions (15 m, 30 m, 70 m, 100 m, 200 m, and 300 m) and found that the largest acceptable DEM resolution was 100 m for a basin with a drainage area of ~300 km². Chaubey et al. [48] found that to get less than 10% bias in Soil and Water Assessment Tool (SWAT) output for total phosphorus (TP), discharge, and NO₃-N simulation, the adequate DEM resolution should be between 100 m and 200 m. While these studies provide valuable information about the effects of the computational cell size and the DEM dataset resolution individually, the joint effect of DEM resolution and computational cell size on modelling performance is yet to be studied.

Finer DEM resolution and finer computational cell size can better capture the basin microtopographic features in hydrological modelling. However, it also means larger computational complexity and longer calculation time for distributed hydrological models. Thus, researchers must find a balance between modelling efficiency/accuracy and the required DEM resolution and the computational cell size. On the premise of ensuring the results of the runoff simulation, the shortest possible runoff simulation time and the smallest workload are pursued. This study, therefore, aims to investigate how DEM resolution and computational cell size jointly affect the performance of routing methods in terms of discharge simulation in distributed hydrological models. Two representative routing methods are selected in this study, including the Liner-reservoir-routing method (LRR), and the improved aggregated network-response function routing method (I-NRF) [49]. LRR is a typical Cell-to-Cell routing method and I-NRF is a typical Source-to-Sink routing method. Two rainy basins with different areas are used in this study. This study can help researchers to choose the most appropriate modelling approach based on the existing data and computation services when predicting drought and flood disasters. In the rest of this paper, Section 2 describes the study area and the data, Section 3 presents a detailed description of the methods and Section 4 presents our results. Then, the conclusions are drawn in Section 5.

2. Study Area and Data

Two basins with different areas are selected to test the generality of the results, namely the Jianxi basin and Shizhenjie basin in southeast China. The Jianxi basin is approximately twice the size of the Shizhenjie basin. The Jianxi River basin (26°31′–28°31′ N, 117°31′–119°00′ W) is the largest tributary in the upper reaches of the Minjiang River, with a drainage area of 14,898 km² above the Qilijie station. Qilijie station is located at the outlet of Jianxi River. The elevation of the Jianxi basin ranges from 166 m to 2008 m (a.m.s.l.). The Lean River has a drainage area of 8294 km² above Shizhenjie station. Shizhenjie station is located downstream of Lean River. In the Shizhenjie basin (28°56′–29°58′ N, 116°94′–118°23′ W), the elevation ranges from 10 m to 1605 m (a.m.s.l.). The drainage net and terrain of our two study areas are shown in Figure 1.

DEM datasets with four different resolutions are used: 1000 m [50], 500 m (http://www.resdc.cn/data.aspx?DATAID=123 (accessed on 14 January 2022)), 250 m (http://www.resdc.cn/data.aspx?DATAID=123 (accessed on 14 January 2022)), and 90 m (http://www.resdc.cn/data.aspx?DATAID=123 (accessed on 14 January 2022)). The daily precipitation, air temperature, and relative humidity data are provided by the National Meteorological Information Centre in China. The daily discharge data of the Qilijie hydrological station in the Jianxi basin from 2000 to 2010 was obtained from the Fujian Hydrology Bureau, and the daily discharge data of the Shizhenjie hydrological station from 1983 to 1993 was provided by the Bureau of Hydrology of the Changjiang Water Resources Commission.

The Long-Range Dependence of the hydroclimatic variables is important in hydrological modelling [51]. Thus, the statistical characteristics of the daily precipitation, air temperature, relative humidity and discharge for the tested period in the study area are calculated and shown in

Table 1. The statistical characteristics of the hydroclimate data in the study area.

		Mean	Standard Deviation	Skewness Coefficient	Kurtosis Coefficient
Jianxi basin	discharge	463 m3/s	753.07	7.88	95.83
	precipitation	1741 mm/year	10.55	4.05	25.20
	air temperature	18 °C	7.66	−0.30	1.94
	relative humidity	78	9.78	−0.22	2.84
Shizhenjie basin	discharge	297 m3/s	515.53	5.40	43.77
	precipitation	1868 mm/year	11.50	3.91	23.62
	air temperature	17 °C	8.81	−0.13	1.79
	relative humidity	79	10.06	−0.31	2.91

. Four important statistical moments [51] are used in this study: the mean, standard deviation, skewness coefficient and kurtosis coefficient. Figure 2 shows the change in the runoff coefficient in the tested period.

3. Methods

3.1. Routing Methods

3.1.1. The Improved Aggregated Network-Response Function Routing Method

I-NRF [49], the improved version of NRF, achieves accurate simulation of discharge results in many basins. It is a typical Source-to-Sink routing method. The I-NRF routing method relies on slope and routing length data extracted from the grids of the DEM dataset:

$$t_i=\left\{\begin{array}{cc}{\displaystyle \sum _{i=1}^k\frac{l_i}{v_{45}\cdot {\left(\tan ({\alpha }_i)\right)}^b}}& {\alpha }_i>{\alpha }_0\\ {\displaystyle \sum _{i=1}^k\frac{l_i}{v_{45}\cdot {\left(\tan ({\alpha }_0)\right)}^b}}& {\alpha }_i\le {\alpha }_0\end{array}\right.$$

(1)

where $t_i$ (minute) is the travel time between each grid and the outlet, k (-) is the total number of grids on the flow path, ${\alpha }_i\left(°\right)$ is the slope, $l_i$ (m) is the routing length between each grid and the downstream grid, and ${\alpha }_0$ ($°$) is the slope threshold that prevents the occurrence of an infinite $t_i$ parameter value. Parameter $v_{45}$ (m/s) is the wave velocity in the grid with a slope of 45$°$ and b (-) is a parameter that reflects how sensitive the wave velocity is to the slope. These two routing parameters $v_{45}$ and b both require special calibration.

3.1.2. Linear Reservoir Routing Method

Based on the characteristic river length theory, the linear reservoir routing (LRR) method, which is a Cell-to-Cell routing method, is popular and reliable [18]. The LRR algorithm is defined as:

$$Q_2=\frac{0.5\Delta t}{0.5\Delta t+K_l}(I_2+I_1)+\frac{-0.5\Delta t+K_l}{0.5\Delta t+K_l}{Q}_1$$

(2)

$$K_l=\frac{l}{v}$$

(3)

where I₁ (m³/s) and Q₁ (m³/s) are the inflow and outflow at the start time, I₂ (m³/s) and Q₂ (m³/s) are the inflow and outflow at the end time, and $\Delta t$ (s) is the calculation time. K_l (s) is a constant that represents the time needed to travel through the characteristic river length, l (m) is the length of the river section, and v (m/s) is the wave velocity of the flood. Of these parameters, only v requires calibration.

The key routing parameters for the I-NRF and LRR routing methods are shown in

Table 2. Descriptions of routing parameters in the I-NRF routing method and LRR routing method.

Routing Method	Parameters	Explanation	Prior Range for Calibration
I-NRF	$v_{45}$ (m/s)	the wave velocity of a grid whose slope is 45°	4, 5, 6, 7, 8, 9, 10	The combination of all values generate 49 routing parameter-value sets
	b (-)	power exponent reflecting how sensitive is the v to slope	0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5
LRR	$v$ (m/s)	The wave velocity	Jianxi: 0.6, 0.8, 1.0, 1.2, 1.4, 1.6 Shizhenjie: 0.2, 0.4, 0.6, 0.8, 1.0, 1.2	There are 6 routing parameter values

.

The LRR routing method needs to satisfy the Courant criterion [52]. The Courant number is defined as:

$$Courant=\frac{u\times dt}{dx}$$

(4)

where Courant is Courant number, u is the wave velocity, dt is the time step, dx is the cell size.

To avoid getting wrong discharge simulation performance, the computational Courant number should satisfy the following criteria:

$$Courant<1$$

(5)

In this study, daily discharge at the hydrological station is used to compare with the simulated e. However, the running time step is set to be 1/20 day when the LRR routing method is used for all computational cell sizes.

The maximum wave velocity and minimum computational cell size that appeared in this article are used to calculate the maximum Courant number in this article.

$$Couran{t}_{max}=\frac{u\ast dt}{dx}=\frac{1.2 \mathrm{m}\times {\mathrm{s}}^{-1}\times 24\times 60\times 60 \mathrm{s}/20}{8180 \mathrm{m}}=0.84$$

(6)

The Courant number is always less than 1 so that the change in the computational cell size does not cause the model to violate the courter number with the LRR routing method.

3.2. Runoff Generation Models

3.2.1. SIMHYD

In this study, two runoff generation models are selected to test the generality of the result. The conceptual daily rainfall-runoff model SIMHYD [53,54] is employed as one of our runoff generation models. This model has been previously applied in numerous study areas, including China [53,55,56,57]. For more information on the model algorithm, refer to Chiew et al. [53]. The details of the seven model parameters in the SIMHYD model are shown in

Table 3. Descriptions of the SIMHYD model parameters.

Parameters	Description	Prior Range for Calibration
INSC (mm)	Interception store capacity	[0.05 10]
COEFF (mm)	Maximum infiltration loss	[0 500]
SQ (-)	Infiltration loss exponent	[0 10]
SMSC (mm)	Soil moisture store capacity	[0 1000]
SUB (-)	Constant of proportionality in interflow equation	[0 1]
CRAK (-)	Constant of proportionality in groundwater recharge equation	[0 1]
K(-)	Base-flow linear recession parameter	[0 1]

.

3.2.2. Water and Snow Balance Modelling System (WASMOD)

The second model used in the study is the water and snow balance modelling system (WASMOD) [58]. WASMOD is a popular model that has been used in many study areas worldwide [59,60,61,62]. The version of WASMOD used in this study is that of Widén Nilsson et al. [62], which consists of only three parameters and no snow module. The parameters for this version of WASMOD are listed in

Table 4. Descriptions of the WASMOD model parameters.

Parameters	About	Prior Range for Calibration
a4 (-)	Actual evaporation	[0.1 0.999]
c1 (1/mm)	Fast runoff	[0 0.1]
c2 (1/mm)	Slow runoff	[0 0.1]
c4 (-)	Potential evaporation	[0.1 0.999]
c5 (-)	Precipitation	[0.5 1.5]

. Compared to the equations listed in Widén Nilsson et al. [62], two additional parameters are added. Parameter $c_5$ corrects the precipitation bias that arises due to a lack of snow measurements that occur occasionally at higher altitude areas and $c_4$ is described in the potential evaporation calculation equation ($c_4$ is not needed when $pe{t}_t$ data are available):

$$pe{t}_t=c_4\cdot {(\max (tm{p}_t+25,0))}^2\cdot (100-r{h}_t)$$

(7)

where $pe{t}_t$ (mm$·$day⁻¹) is the potential evaporation, $tm{p}_t$ ($℃$) is the temperature, and $r{h}_t$ (-) is the relative humidity, $c_4$ (mm day⁻¹ $℃$⁻²) is the parameter.

3.3. Different Computational Cell Sizes and Interpolation Method

Following the boundary of the study basin, computational cells for routing methods are generated. The computational cell sizes range from 5 arc-min to 60 arc-min with an interval of 1 arc-min, resulting in a total of 56 cell sizes. For illustration purposes, the computational cells of Jianxi basin at the resolution of 10′ and 30′ are shown in Figure 3. The computational cell is constructed according to the following principle: the leftmost longitude of the cells is the same as the integer part of the leftmost latitude of the basin boundary and the bottom latitude of the cell is the same as the integer part of the bottom latitude of the watershed boundary. To maintain the water balance, the area of each cell is the area of the basin belonging to this cell. It can avoid the area errors for the boundary cells [18]. After generating the computational cells, extract the flow net at the same resolution of computational cell size using the network scaling algorithm (NSA) [63] (Figure 3).

In order to match the resolution of the meteorological data recorded at the hydrological station with the various computational cell sizes, the Inverse Distance Weighted (IDW) method [64] is employed to interpolate the data to the grid with a resolution of 1 km. For each grid, the distance from the grid center to all hydroclimatic stations is calculated. The four stations closest to the center of the grid are selected to calculate the hydroclimatic data values for this grid. The equations are as follows:

$$weight(i)=(1/dis(i))/{\displaystyle \sum _{i=1}^4(1/dis(i))}$$

(8)

$$data=data\_observation(i)\times weight(i)$$

(9)

where i is the corresponding number of the 4 stations which are the closest stations to the grid center, weight is the weight of the station data, data_observation is the measured data of the station, and data_g is the hydroclimatic data of the grid.

Using the linear average method, the interpolated meteorological data are then aggregated into the gridded data with the resolution that matches that of the computational cell size. For each cell, the hydroclimatic data value is calculated based on the grids which are belonging to the cell.

$$data\_c=\frac{{\displaystyle \sum _{j=1}^{n}data\_g(j)}}{n}$$

(10)

where data_c is hydroclimatic data of the cell, j is the corresponding number of the grids belonging to the cell.

3.4. Calibration Algorithms

Two calibration algorithms are used in this study: the Monte Carlo algorithm [65] and the Covariance Matrix Adaption Evolution Strategy (CMA-ES) [66]. Since the available data are not long enough, the calibration period length is set as 7 years and the validation period is set as 4 years according to Klemeš [67]. Besides, the warming-up period equals 1 year in this study. The detailed information on the runoff generation and routing parameters that require calibration are listed in Table 2, Table 3 and Table 4.

The Nash–Sutcliffe efficiency (NSE) [68] is used for model calibration criterion, percentage bias (PBIAS) [69] is also used for results evaluation and comparison.

$$NSE=1-\frac{{\displaystyle {\sum }_{i=1}^n{(O_i-S_i)}^2}}{{\displaystyle {\sum }_{i=1}^n{(O_i-\overline{O_i})}^2}}$$

(11)

$$PBIAS=\frac{{\displaystyle \sum _{i=1}^n(O_i-S_i)}}{{\displaystyle \sum _{i=1}^n{O}_i}}\times 100$$

(12)

where $O_i$ (m³/s) is the observed discharge, $\overline{O_i}$ (m³/s) is the mean observed discharge, $S_i$ (m³/s) is the simulated discharge, and n (day) is the length of the time step.

3.5. Experimental Design

There are a total of 1792 experiments in this study based on the cross combination of two study basins, two routing methods, two runoff generation models, four DEM resolutions and fifty-six computational cell sizes.

3.5.1. Optimal Routing Parameters

In this study, according to the characteristics of the routing method, the Monte Carlo algorithm is used to calibrate the optimal I-NRF routing parameters, and the CMA-ES algorithm is used to calibrate the optimal LRR routing parameters.

The routing parameters are mostly related to basin physical characteristics such as the topography and the underlying surface. As such, it is assumed that the I-NRF and LRR routing parameters are independent of the parameters from the runoff generation models.

Calibration of I-NRF Routing Parameters

The calibration process for the I-NRF routing parameters consists of four steps:

• Several initial routing parameter sets are drafted. For the I-NRF routing method, there are two parameters to be calibrated. Table 2 defines 49 initial routing parameter set values. The number of the initial values is m (m = 1~49).
• The optimal simulation performance for each initial routing parameter for a given DEM resolution and a given computational cell size is calculated using Monte Carlo calibration algorithm. The optimal model performance for each parameter set is denoted as $NS{E}_{i,j}^m$, where i corresponds to the DEM resolution and j corresponds to the computational cell size. Then 48 $NS{E}_{i,j}^m$ values using the four different DEM dataset resolutions (i = 1–4) and the 12 computational cell sizes (5 arc-min to 60 arc-min with an interval of 5 arc-min) (j = 1–12) are obtained. In the calibration process, 300 runoff generation parameter arrays are produced as the model input by Latin-Hypercube sampling [70] according to the initial parameter range (Table 4). The applied marginal distribution is the uniform distribution.
• The comprehensive simulation performance for each initial routing parameter set of each DEM resolution ($NS{E}_i^m$) is calculated, where $NS{E}_i^m=({\displaystyle \sum _{j=1}^{12}NS{E}_{i,j}^m})/12$.
• The optimal routing parameter set is selected from the initial routing parameter sets. For each DEM resolution, the optimal routing parameter set for a given basin is selected using $\underset{m=1}{\overset{49}{\max }}(NS{E}_i^m)$ as a criterion and listed in

  Table 5. The optimal routing parameters of 4 kinds of DEM resolutions.

  DEM Resolution (m × m)	Jianxi Basin			Shizhenjie Basin
  	v45 (I-NRF)	b (I-NRF)	v (LRR)	v45 (I-NRF)	b (I-NRF)	v (LRR)
  90 × 90	10	0.35	1.2	4	0.25	0.6
  250 × 250	5	0.2	1.2	6	0.3	0.6
  500 × 500	9	0.3	1.2	8	0.35	0.6
  1000 × 1000	7	0.2	1.2	9	0.35	0.6

  .

Calibration of LRR Routing Parameters

The CMA-ES algorithm is used to identify the optimal routing parameters for the LRR routing method.

The calculation of the optimal LRR routing parameters happens in four stages:

• Several initial routing parameters are drafted. For the LRR routing method, there is only one parameter to be calibrated. Table 2 defines 6 initial values. The number of the initial values is n (n = 1–6).
• The optimal simulation performance for each routing parameter for a given DEM resolution and a given computational cell size is calculated using CMA-ES algorithm, denoted as $NS{E}_{i,j}^n$. With 4 DEM resolutions (i = 1–4) and 12 computational cell sizes (5 arc-min to 60 arc-min with an interval of 5 arc-min) (j = 1–12), 48 $NS{E}_{i,j}^n$ values are obtained. In the calibration process, the initial range of runoff generation parameters is given in Table 4.
• The comprehensive simulation performance ($NS{E}_i^n$) for each parameter of each DEM resolution is calculated. Here, $NS{E}_i^n=({\displaystyle \sum _{j=1}^{12}NS{E}_{i,j}^n})/12$.
• The optimal routing parameter is selected from the initial routing parameters. For each DEM resolution, the optimal routing parameter set for the basin is selected using $\underset{n=1}{\overset{6}{\max }}(NS{E}_i^n)$ as a criterion and the results are listed in Table 5.

3.5.2. Optimal Runoff Generation Parameters

As mentioned above, there are a total of 1792 experiments. Once the optimal routing parameters are obtained and fixed, the runoff generation parameters are calibrated in each experiment using the CMA-ES algorithm [66]. Then the corresponding NSE values are ultimately obtained (Figure 4, Figure 5 and Figure 6) for 1792 experimental simulations.

4. Results

4.1. The Effect of Computational Cell Size on the Discharge Simulation Performance

The NSE values for the various DEM resolutions (x-axis) and computational cell sizes (y-axis) for both routing methods in Jianxi basin and Shizhenjie basin are shown in Figure 4 and Figure 5. To improve readability, line charts of NSE are shown in Figure 6.

Figure 4, Figure 5 and Figure 6 show that the discharge simulation performance curve for the Jianxi basin, as computational cell size changes, is mostly a straight line. While the WASMOD model generates some noise in the NSE values of the Shizhenjie basin (Figure 5a,b and Figure 6b), the NSE values are mostly constant, despite the changing computational cell sizes. In general, the discharge simulation performance with the I-NRF routing method stays stable, regardless of the computational cell size. Gong et al. [18] discussed in detail why NRF discharge simulation performance is scale-independent. I-NRF, the improved version of NRF, retains the scale-independency characteristic of NRF. This scale-independency is produced by the algorithm, in which the discharge simulation performance is obtained by the convolution of the runoff generation result and the cell response function (CRF). The CRF is aggregated from the pixel response function (PRF), which is determined by the slope and the routing length extracted from the DEM grid when the wave velocity is fixed. Besides, there is no loss of information in the process of converting the PRFs to CRFs. Thus, the CRFs contain all the information from the PRFs and the DEM. As long as the DEM resolution stays the same, this algorithm generates stable results that are insensitive to the changes of computational cell size.

The performance of the LRR routing method is, however, different from that of the I-NRF routing method, as it fluctuates as the computational cell size changes (Figure 4c,d Figure 5c,d and Figure 6c,d). The PBAIS shows the same trend as NSE. The computational cell size affects the LRR PBAIS more than I-NRF PBIAS. Furthermore, the fluctuations observed in the NSE values in the Jianxi and Shizhenjie basins are not identical. For Jianxi basin, the NSE curve fluctuates greatly, from 20′ to 21′ and from 30′ to 31′, etc., for every DEM resolution (Figure 4 and Figure 6). For the Shizhenjie basin, the NSE curve fluctuates from 18′ to 19′ and from 28′ to 29′, etc., with DEM of 90 m. It fluctuates greatly from 19′ to 20′ and from 28′ to 29′, etc., with DEM of 250 m, 500 m and 1000 m (Figure 5 and Figure 6).

The sensitivity of the LRR routing method to the computational cell size arises because the flow net changes with the computational cell size [18]. In fact, the flow net experiences distinct structural changes as the computational cell size changes. For example, the NSE curve fluctuates greatly, from 20′ to 21′ and from 30′ to 31′, etc., for every DEM resolution in Jianxi basin (Figure 4 and Figure 6). The differences of computational cell size between 20′ and 21′ or 30′ and 31′ are very small, but the resulting NSE differences are significant. To better understand these turning points, the extracted river networks for computational cell sizes of 20′, 21′, 30′, and 31′ in the Jianxi basin are plotted in Figure 7. From 20′ to 21′ and from 30′ to 31′, minor changes in the computational cell size lead to a different allocation of the outlet grid cells, which ultimately results in a significant difference in the flow net. Besides, the LRR routing method is based on the flow net at a given computational cell size. As such, a small change in the computational cell size leads to a large difference in the NSE value.

Additionally, the LRR NSE values are also sensitive to the computational cell size because smaller computational cell sizes correspond to a finer and different flow net (Figure 3). The high-resolution flow net not only better captures the shape of the real river network, but also better simulates the length of the routing path. In the LRR method, the length of the routing path affects the deformation and the smoothing of the flood. A coarser computational cell size reduces the flow net between cells into a straight line, which explains why the NSE value is higher and PBIAS is smaller when the computational cell size is small. The result that the LRR routing method is sensitive to the computational cell size confirmed the result from Arora et al. [45]. In that study, a finer computational cell size resulted in better performance.

Whether the NSE of the LRR and I-NRF routing method is sensitive to computational cell size is not related to the runoff generation model and the specific study basin in this study, indicating that this pattern is universal. Therefore, to accurately simulate the discharge, larger computational cell sizes can be used to save computing time when the I-NRF routing method is used. For the LRR routing method, however, the choice of computational cell size must be an informed decision.

4.2. The Effect of the DEM Resolution on the Discharge Simulation Performance

Box plots of NSE value variations with the computational cell size and the DEM resolution are shown in Figure 8. The impact of the DEM resolution on the discharge simulation performance is slightly larger for the I-NRF routing method than it is for the LRR routing method. For the I-NRF routing method, the use of different DEM resolutions leads to some differences in modelling performance, for both runoff generation models and for both study basins; from these observations, it can be inferred that the DEM resolution directly affects the NSE values, to some extent (Figure 8c,d,g,h). For the LRR routing method, the average NSE for all computational cell sizes is similar when DEM resolution changes, regardless of the runoff generation model or the study basin. As such, it can be concluded that the DEM resolution has little effect on the overall LRR discharge simulation performance (Figure 8a,b,e,f). The laws analyzed from changes in PBIAS are the same as those analyzed from changes in NSE.

With the I-NRF routing method, different DEM resolutions produce variations in the overall discharge simulation performance (Figure 6 and Figure 8). However, in this study, the difference in the NSE values generated by the different DEM resolutions is less than 0.03 (Figure 6), while the difference in the NSE values due to the implementation of the two runoff generation models is greater than 0.03 (Figure 6). Besides, as shown in Figure 8, there is no consistent correlation between the DEM resolution and the overall discharge simulation performance for the I-NRF routing method. In the I-NRF routing method, an accurate simulation of discharge value requires an accurate simulation of the routing time. As per Equation (1), the routing time $t_i$ is calculated using the slope and the routing length. The slope and routing length both increase with a finer DEM resolution [46]. In practice, it is difficult to predict whether the $t_i$ value is closer to or further from the true routing time value when the DEM resolution becomes finer (Equation (1)). Consequently, a finer DEM resolution does not necessarily correspond to higher NSE values, a fact that was observed in our study. In the Jianxi basin, the worst discharge simulation is generated with a DEM resolution of 1000 m $\times$ 1000 m; in the Shizhenjie basin, the worst discharge simulation is generated with a DEM resolution of 90 m × 90 m. Figure 3, Figure 4 and Figure 6 also show that DEM resolution and computational cell size have almost no joint influence on the result of I-NRF.

For the LRR routing method, the NSE difference introduced by DEM resolution is much smaller than that of computational cell size (Figure 6 and Figure 8a,b,e,f). Finer DEM resolution does not result in better model performance (Figure 8a,b,e,f). As discussed previously, the LRR routing method is based on the flow net at a given computational cell size. In addition, even the coarsest DEM resolution (1000 m) is still much smaller than the smallest computational cell size (about 8180 m). Thus, the DEM slope and routing length are largely irrelevant. Combining Figure 8 with Figure 6, it can be seen that the NSE with a DEM resolution of 500 m is different from others at the computational cell size of 26′ and 39′ in the Jianxi basin. Furthermore, for the Shizhenjie watershed, the NSE curves corresponding to each DEM resolution do not overlap. To understand the difference caused by the DEM resolution, the flow nets at the computational cell size of 39′ with DEM resolutions of 1000 m, 500 m, 250 m, 90 m in the Jianxi basin are drawn in Figure 9. It can be seen that the flow net with the DEM resolution of 500 m is different from the others. That is the reason why NSE with the DEM resolution of 500 m is different from others at the computational cell sizes of 39′. Now it is known that the DEM resolution and the computational cell size affect the upscaled flow net at the same time. On the one hand, when the DEM resolution is different, the flow direction identification of some computational cells may be different, which will cause disturbance to the basin flow net. On the other hand, the computational cell size will also affect the basin flow net, as described in Section 4.1. Because the flow net directly affects the discharge simulation performance in the LRR routing method, it can be proposed that the DEM resolution and the computational cell size jointly affect the simulated discharge when the LRR routing method is implemented. However, the joint effect is very small and can be ignored in the actual simulation (Figure 6 and Figure 8).

In this paper, both coarser and finer DEMs are sufficient for hydrological modelling with those conceptual models. However, such conclusions are not suitable for most hydraulic models. When using a hydrodynamic model, caution is still required regarding the choice of DEM resolution. For example, Horrit et al. [71] studied the effect of DEM resolution on the LISFLOOD-FP model, using a 1D description of channel flow. Models with resolutions varying from 10 to 1000 m are compared. The results showed that the model achieves the best simulation of the inundated area when using a 100 m DEM. Cook and Merwade [72] found that using a 1D HEC-RAS model, and a 2D FESWMS model, the flood inundation area reduces with improved horizontal resolution and vertical accuracy in the topographic model.

In order to facilitate researchers to choose the appropriate routing method when conducting hydrological simulations, this paper lists the time duration of the two models. Taking the case of using the WASMOD runoff generation model and Jianxi basin as an example, the running times of the two routing methods are compared (

Table 6. The time duration of the hydrological model with I-NRF or LRR routing method.

		I-NRF	LRR
Jianxi Bain (1000 m)	5 arc-min	1768 s	2450 s
	60 arc-min	165 s	225 s

). According to the above discussion, for both routing methods, a finer DEM cannot guarantee a better simulation effect. Therefore, the time duration of the model under different DEM resolutions is not discussed here. It should be noted that the time duration of the models listed in this article does not include pre-processing and preparation of data, but includes parameter optimization of the model. The computer used to count the running time is as follows. OS: Windows 10 Pro 64-bit; CPU: Intel Core i7-9700 @ 3.00GHz eight-core; Memory: 16 GB; Sofeware: MATLAB R2020b.

It can be seen that the time duration of both routing methods is within the acceptable range, and the time duration of I-NRF is slightly shorter. However, a finer computational cell size leads to a huge increase in time duration. Thus, a finer computational cell size means a longer time duration and better simulation performance for the LRR routing method, while coarser computational cell size can be used with the I-NRF routing method to save time. Therefore, this paper proposes to choose the I-NRF routing model when conducting watershed-scale hydrological simulations.

4.3. Analysis of Discharge Duration Curves

To further analyze the effects of the DEM resolution and the computational cell size on the simulated discharge for the two routing methods, the corresponding discharge duration curves are calculated for each DEM resolution and computational cell size. The discharge duration curve is a statistical curve that represents the proportion of days in which the flow exceeds a specified threshold value in a certain time period. Figure 10 and Figure 11 plot the observed discharge duration curve and the simulated discharge duration curve clusters for the various DEM resolutions and computational cell sizes. The colors of the simulated discharge duration curve clusters correspond to the different DEM resolutions. The ordinate is the daily flow and the abscissa is the cumulative proportion. To more easily visualize the instances of high flow predicted by the hydrological models, a horizontal logarithmic axis is applied to this plot. Smaller values on the horizontal axis correspond to larger flows.

As shown in Figure 10 and Figure 11, both DEM resolutions and computational cell size influence the simulated discharge. The effects of the computational cell size and the DEM resolution on the simulated discharge accuracy are mainly reflected in the high-flow range (i.e., the first 1% of the discharge). As shown by the discharge duration curves of the two runoff generation methods, it can be seen that the WASMOD model tends to underestimate the high flow, and SIMHYD can simulate high flow better than WASMOD in the study area.

4.4. Analysis of Runoff Generation Parameters

The optimal runoff generation parameters in the Jianxi and Shizhenjie basins, identified from a combination of 4 DEM resolutions and 56 computational cell sizes, are shown in Figure 12 and Figure 13.

In the Jianxi basin, the optimal runoff generation parameters with the I-NRF routing method vary with the different DEM resolutions (Figure 12). Furthermore, when the runoff generation parameters with the I-NRF routing method change with the computational cell size, parameters $a_4$, $c_1$, and $c_2$ in the WASMOD model are more stable than the parameters in the SIMHYD model, a difference that may be caused by the equifinality effect that arises due to the higher number of parameters in the SIMHYD model. The equifinality of model parameter sets is a common phenomenon. There are records of this phenomenon in Beven [73,74]. The optimal runoff generation parameters with the LRR routing method appear to be largely insensitive to the DEM resolution. With the LRR routing method, the turning points in the curves in Figure 12 correspond somewhat to those in Figure 6c. Additionally, they show the same or the opposite trend to the NSE curves. Some parameters are sensitive to the computational cell size, while others are not. In the Jianxi basin, parameters SUSC, K, $c_1$, and $c_2$ exhibit opposing NSE trends. While, parameters SQ, SMSC, SUB, and K in the SIMHYD model and parameters $a_4$, $c_1$, and $c_2$ in the WASMOD model, parameters $a_4$, SUB, and SQ, exhibit similar trends with NSE (Figure 12 and Figure 6c).

The results for the Shizhenjie basin are shown in Figure 13. Just as they did in the Jianxi basin, the optimal runoff generation parameters with the I-NRF routing method vary with the DEM resolutions. The runoff generation parameters with the I-NRF routing method vary much more widely with the computational cell size in the Shizhenjie basin than they do in the Jianxi basin. Besides, the runoff generation parameters with the LRR routing method vary much more widely with the DEM resolutions in the Shizhenjie basin than they do in the Jianxi basin. With the LRR routing method, the parameters $c_1$ and SUB follow the same NSE trends, while parameters K and $c_2$ exhibit opposing NSE trends as computational cell size changes in the Shizhenjie basin (Figure 13 and Figure 6d).

The uncertainty in simulation performance includes the input forcing, boundary condition data, model structure and so on [75,76]. It will be necessary to study how they affect the overall uncertainty in the on-going research.

5. Conclusions

Previous studies have reported how the DEM resolution and computational cell size affect the modelling results separately. However, too little work has been devoted to how they jointly affect the routing method performance. This study sought to understand how the DEM resolution and the computational cell size jointly affect the simulated discharge results produced by the I-NRF and LRR routing methods. Two runoff-generation models were used to test the generality of the result. In these models, computational cell sizes that range from 5 arc-min to 60 arc-min, with an interval of 1 arc-min and DEM resolutions of 90 m, 250 m, 500 m, and 1000 m are used. The conclusions are as follows:

1.

The DEM resolution has a larger impact on the I-NRF discharge simulation performance than it does for the LRR method. However, the overall I-NRF discharge simulation performance is not consistently correlated with DEM resolution. As the computational cell size increases, the I-NRF discharge simulation performance is relatively stable because the I-NRF method is largely independent of the computational cell size. The DEM resolution and the computational cell size exert no joint influence on the I-NRF discharge simulation performance.

2.

The DEM resolution has little effect on the LRR discharge simulation performance. The LRR discharge simulation performance changes drastically with the computational cell size for both runoff generation models. In general, finer computational cell size leads to better results, while the LRR discharge simulation performance oscillates heavily as the computational cell size increases. The joint effect of DEM resolution and the computational cell size on the performance of the LRR routing method can be ignored in hydrological modelling.

3.

The effects of the computational cell size and the DEM resolution on the simulated discharge accuracy are mainly reflected in the high-flow range. In this study, the performance of SIMHYD is better than that of WASMOD in the high-flow range, when they are combined with both routing methods.

These results will provide guidance as to how hydrological researchers can use their knowledge of the study area and their computational limitations to make reasonable routing method, DEM resolution, and computational cell size choices.