Dynamic Manning's roughness coefficients for hydrological modelling in basins

Most recently used routing models have been based on the Saint Venant equations or their approximations, such as a kinematic wave, noninertia wave, linear diffusion-wave (Wang et al. 2014) and quasi-steady dynamic wave (Reggiani et al. 2014). Manning's roughness coefficient, n, is one of the most important parameters in hydrological calculations representing the loss of energy in open channels. It is commonly used to calculate discharge and flood water elevations (Coon 1995). Usually, the n value is a parameter of the routing module in a hydrological model, Manning's formula is used to calculate discharge and flood levels in equilibrium conditions. The value of n has an important effect on the accuracy of the simulated streamflow. However, n is difficult to obtain in natural basin (channels and slopes) because it incorporates many factors including vegetation density, riverbed irregularity, surface water width and soil component differences, which contribute to the resistance of flow. The values of Manning's roughness coefficients are different in channel and slopes, the main cause being that the landscape surface in slopes is rougher than channel (Cowan 1956; Azamathulla & Jarrett 2013). Hence, substantial attention has been focused on the roughness coefficient calculations in natural basins.

Equation (1) can be used to calculate uniform flow in which the water-surface profile and energy gradient are parallel to the riverbed, and the river cross-section area, hydraulic radius and depth remain constant throughout the river reach (Jarrett 1985).

Considering uncertainties, Cowan (1956) proposed a formula using five parameters (n = (n_b + n₁ + n₂ + n₃ + n₄)·m) to represent the influence of different variable factors, where n_b is a base value of n for a straight, uniform, smooth channel in natural materials; n₁ is a correction factor for the effect of surface irregularities; n₂ corrects for variations in the shape and size of the channel cross section; n₃ corrects for obstructions; n₄ corrects for vegetation and flow conditions; and m is a correction factor for meandering. Petryk & Bosmajian (1975) developed a method to derive n based on the vegetation density for a densely vegetated flood plain. The vegetation density has a strong impact on the roughness of a channel (Li et al. 2014). An equation was derived that incorporates the influence of the stem density on the flow resistance (Noarayanan et al. 2012); it was found that the energy loss of the head due to friction is caused by both vegetation and the side walls. Water flows through different river cross sections with different hydraulic radiuses. When the water depth is far less than the water flow width, the hydraulic radius is approximately equal to the water depth. Thus, n varies dynamically with the streamflow and changes in water depth. The value of n decreases with increasing depth of flow and varies directly with slope for high-gradient streams (Azamathulla & Jarrett 2013). However, the water depth is difficult to measure for all cross sections in an entire basin. Some equations for calculating n are based on an analytical model and solved by iteration (Li & Zhang 2001). The characteristic size of streambed particles has an influence on the roughness. Some studies investigated the relationships between the distribution of the particle sizes and n (Limerinos 1970).

The roughness coefficient is not a constant because the water depth and vegetation density change during different seasons in a basin. However, current hydrological models often apply a static n in calculations, without considering its dynamic changes. The roughness coefficient is a sensitive parameter in hydrological models (Ye et al. 2013). By contrast, some studies considered it to be an empirical parameter and calibrated the hydrological model by modifying the n value until the simulated discharge reasonably matched the observation (Mahmoudi et al. 2015). The n value based on more physical mechanisms can decrease the model uncertainty and parameter sensitivity (Aronica et al. 1998). There are obvious errors in the simulation for sub-basins, although the discharge simulation is satisfactory in large basin outlets because of equifinality for different parameter sets (Beven 2006). To simulate a flash flood or low flow and avoid parameter over-optimization in each sub-basin, dynamic roughness coefficient estimation (high-accuracy n) is absolutely necessary even for a large scale basin.

To manage the n value over a relatively accurate range, some studies have divided river cross sections into different strips with uniform roughness coefficients based on the land cover classification, which were from aerial photos and field surveys, to calibrate the n value by comparison with the observed water levels (Kovacs et al. 2006; Tóth 2009; Gichamo et al. 2012).

In this study, we tried to propose dynamic Manning's roughness coefficients (DMRC) and improve the routing of a hydrological model in basins.

The structure of the paper is as follows: immediately below the hydrological model and statistical method are described; a section introducing the data and study domain follows; next, a section presents the results and discussion; and the final section provides the conclusions.

A new equation is proposed to simulate n using statistical analysis of the estimated n values from Manning's formula (Equation (1)) combined with the LAI, soil components, water depth and water surface width of a river. The new equation was then applied to improve the routing module of the Distributed Time-Variant Gain Hydrological Model (DTVGM) (Xia et al. 2005) (Figure 1). The DTVGM inputs are the precipitation, temperature and types of vegetation. The DTVGM outputs are discharge and evapotranspiration in sub-basins. The DMRC model is a module of DTVGM. The DMRC output is the roughness coefficient in the sub-basins. The DMRC inputs are the LAI, soil structure and water flow cross-section area (m²), which is from the routing model output. The modelling time step occurs daily.

The DTVGM is a daily distributed hydrological model that incorporates runoff and routing modules. Coupling the advantages of both nonlinear and distributed hydrological models, the DTVGM can simulate variable hydrological processes under complex environmental conditions. It has been applied to two different cases: the Heihe River basin, which is in the arid and semi-arid region of north-western China; and the Chaobai River basin, which is in the semi-humid region of northern China (Xia et al. 2005). Multiple versions of the DTVGM have been proposed, and the model still undergoes continuous improvement and innovation.

The n coefficient is often a constant in traditional hydrological models and land surface models. However, the most optimal n for a model is not always equal to a constant, it changes dynamically with the changing water flow cross-section area and vegetation in improved hydrological models.

In prior research, the n value in Manning's formula was a function of the vegetation density, soil components, water depth and surface width of rivers (Moharana & Khatua 2014). In this study, we analysed the relationships between n and the water depth and surface width using linear regression. An optimal equation for calculating n was then derived based on multiple nonlinear regression.

Multiple regression analysis is an appropriate method when the research problem includes one unique metric-dependent variable and a dependent-variable that is related to more than one metric-independent variable (Hair et al. 2006). The general purpose of multiple regression analysis is to learn about the relationship between several independent or predictive variables and a dependent or criterion variable (Enayatollahi et al. 2014). Additionally, an advantage of multiple nonlinear regression is that a very high order multivariate should be able to approximate complex multivariate functions (Cogger 2010).

Experimental results have indicated that Manning's roughness coefficient n increases with the increasing vegetation density, which leads to higher flow resistance in natural channels and hillsides (Li et al. 2014). The vegetation density is difficult to obtain in basins. However, the leaf area index (LAI) can be measured by remote sensing. Instead of the vegetation density, the LAI was used in this study to develop a formula for calculating the n. In the simulation process, the LAI was set to (LAI+ 1) to avoid calculation errors when the LAI equals zero in bare areas.

Stream bed roughness is associated with bed material, especially the particle size and distribution (Limerinos 1970). We selected a soil particle (clay, loam and sand content) dataset for the bed material information.

Equation (13) was used in the routing process in an entire basin; n was separately calculated in each sub-basin during each time period. The water flow cross-sectional area can be computed using Newton iterations (Equation (12)) in both channels and slopes. A continuous and consistent LAI (Yuan et al. 2011) is composited every 8 days at 1-km resolution. The LAI was interpolated into the sub-basins. The LAI in the channels and slopes are equal to the LAI that locates in sub-basins. The A and LAI vary across time and space.

The p₁, p₂ and p₃ parameters were optimized by 1stOpt (http://www.7d-soft.com/). 1stOpt is a robust, easy-to-use and powerful optimization tool that is widely used in various engineering fields. Based on the optimization software package 1stOpt, multiple nonlinear regression can be easily established and solved (Wang et al. 2009; Hu et al. 2011). The 1stOpt search of optimization capability is stronger than that of other simulation software because it can find relatively accurate results from any initial value (Tang et al. 2008). The Levenberg–Marquardt + Universal Global Optimization of 1stOpt was applied to optimize a set of parameters in Equation (2).

There are three steps to calculate dynamic n: (1) reading the soil structure and LAI at t time, (2) deriving the water flow cross-sectional area using Newton iterations (Equation (13)) based on n_t,k−1, (3) using Equation (2) to calculate n_t,k for the next iteration. Therefore, n is dynamic in the improved routing module.

rBias refers to the correspondence between the average simulated value and the average observations. R checks the fitness of the simulated and observed values. The RMSE is used to measure the deviation between the observation and simulated values. If NSE is close to 1, the model is considered to have good quality and be highly reliable. If it is close to 0, the simulation is close to the observed value of an average level and the overall results are legitimate, but a bad simulation error has occurred. If the NSE is smaller than 0, the model is not reliable.

The Yellow River is the second-longest river (Figure 2) in China, originating from the Tibetan plateau. It passes through the northern semi-arid region, crosses the Loess Plateau and finally discharges into the Bohai Gulf (Yang et al. 2004). The main course of the Yellow River flows for approximately 5,464 km and has a drainage basin area of approximately 752,443 km². The Yellow River was divided into 3,316 sub-basins, and each sub-basin area is greater than 100 km². The climate conditions vary from cold to temperate zones and change from arid and semi-arid to semi-humid regions (Cheng 1996). The annual precipitation for the Yellow River is approximately 450 mm. However, the annual evaporation is 1,100 mm.

There are 35 stations with estimated data to simulate the roughness coefficient (Table 1), which include the water surface width (w), discharge (Q), water depth (H), friction slope (s) and roughness coefficients (n) for the period from 2008 to 2012, obtained from hydrological year books of the Yellow River basin. Manning's roughness coefficient is estimated using Manning's formula.

The observed discharge in three hydrological stations was selected to compare the raw (the model with static n) and improved (the model with dynamic n) model results. We selected the three stations because it is very difficult to obtain more observed discharge data in other small catchments. Tangnaihai station has a 121,927 km² drainage area in the upper Yellow River basin; Jingchuan and Yuanjiaan stations have 3,145 km² and 1,661 km² drainage areas, respectively, in the midstream Yellow River basin. Additionally, the observed daily discharge data (from 1998 to 1999) of Tangnaihai station, Yuanjiaan station and Jingchuan station were collected from the Yellow River Conservancy Commission (YRCC).

Daily precipitation and potential evaporation data from 1955 to 2012 were provided by the China Meteorological Administration. There are 983 precipitation gauges and 839 evaporation gauges in continental China. A continuous and consistent LAI (Yuan et al. 2011) is composited every 8 days at 1-km resolution, and a conterminous China multi-layer soil particle-distribution dataset (clay, loam and sand content) was developed by Shangguan et al. (2012) with a resolution of 1 km × 1 km.

The daily precipitation and evaporation data, LAI (Figure 3) and Yellow River soil particle-distribution (Figure 4) were interpolated into the 3,316 sub-basins. Figure 3 shows that the LAI on the slope is much larger than that from the main streams. The LAI value is much larger in the upstream and downstream in July. The main causes are that approximately 70% of the total annual rainfall is restricted to the summer (June–September) (Ye et al. 2006b), and there is a desert in the middle stream in the Yellow River basin.

From the definition of the roughness coefficient, we know that the roughness coefficients were determined by the soil, specific discharge conditions, and size, shape and types of vegetation that line the bed and sides of the channel. Therefore, we selected the composition of the soil (sand, loam and clay), LAI and water flow cross-section area to estimate the n value. We collected 1,177 estimated roughness coefficients and corresponding water depths, surface widths, LAIs and soil data from 35 hydrostations. However, the samples were not sufficient to quantify the effect of soil and vegetation in 35 stations. Li et al. (2014) and Noarayanan et al. (2012) studied vegetation–flow interactions under laboratory conditions. Their study showed that Manning's vegetation roughness coefficient due to vegetation resistance increased with the increasing vegetation density. We used their results in our study. The roughness coefficient is closely related to the particle size and distribution (Limerinos 1970). Limerinos (1970) reported that the large particle size increased the roughness coefficient by using 11 catchment data sites in the USA.

To explore the relationship between the n value and water depth d and surface width w, scatter diagrams were produced to support the statistical analysis.

The observed d and estimated n values for six stations are shown in Figure 5. Obviously, there is a negative relationship between d and n for the Tongren, Huangyuan, Xining, Ledu and Dongjiazhuang stations, meaning that Manning's n decreases with an increase of water depth. This means that less energy was consumed as the flow depth increased in open channels. This behaviour was also confirmed by Fippin-Dudley et al. (1998). However, this negative relationship is not shown for Damitan station. The main cause is that a decreasing n value is a comprehensive function of multiple factors, such as the riverbed irregularity, vegetation density and particle sizes in a riverbed. The riverbed has been cleared and fixed at Damitan station such that the roughness is more like constant at Damitan station.

Figure 6 shows the variation of n with the water surface width. The n value decreases with increasing water surface width, meaning that broader water surface width will generate lower flow resistance.

The hydrological stations shown in Figure 2 are labelled, starting from upstream of the Yellow River to the downstream. There are only 1,177 estimated n values by Manning's formula in 35 stations. The data period is from 2008 to 2012. The period is too short, and we wanted to verify the improved method in different locations. Therefore, the observed data of 21 stations (A1–A21) were used to calibrate the parameters in Equation (2), and the data from the remaining 14 stations (A22–A35) were used to validate the p₁–p₃ parameters.

There are different percentages of clay, loam and sand in each type of soil. With a large particle size, n is large in sand. Furthermore, a high vegetation density contributes to the resistance of flow.

In the DTVGM, each sub-basin is defined as a hydrological unit to calculate the runoff, roughness coefficients and routing.

Table 2 shows the DTVGM parameter values. These parameters were the optimized results from using the observed streamflow. A manual calibration method was used to calibrate the model parameters because the distributed hydrological model takes a long time to run, and automatic calibration would be too time-consuming. During manual calibration, the model was run a few times to ensure that NSE, R and rBias were good. Due to the difference of vegetation, soil and climate between upstream and downstream in the Yellow River basin, the parameters of the runoff module are different at different hydro-stations.

Except for the roughness coefficients, parameters in the raw (static n) DTVGM and improved (dynamic n) DTVGM are the same at the same station. Using the raw and improved DTVGM, we produced raw and improved simulated discharges.

Figures 8–10 show the daily hydrograph of Tangnaihai, Jingchuan and Yuanjiaan stations in 1998 and 1999. In Figure 8(a), the two lines of the raw discharge and improved discharge are almost the same and the difference of error is not obvious, mainly because Tangnaihai station is located in the main channel of the Yellow River, and the catchment area is very large (121,972 km²). Manning's n is greater than true value in main river channels and less than true value in the slopes in the static Manning's n model. The velocity of the flow would decrease if n increases and other conditions do not change, and vice versa. Manning's n is not equal to true value in slopes and river channels, the simulated velocity of the flow is smaller than true value in slopes and is greater than true value in river channels. The total simulated streamflows are similar because the errors are offset in slopes and river channels. Figure 8(b) depicts the absolute error between the observed and simulated discharge. However, the absolute error in Figures 9(b) and 10(b) is more obviously different than that in Figure 8(b). The drainage areas of the Jingchuan and Yuanjiaan catchments are small. The simulated flood peak in the improved model is less than that in the raw model and is more close to the observed discharge in Figures 9(a) and 10(a). The improved model affords a reasonable roughness to the slope, making the new simulated discharge more accurate than the raw simulated discharge in a small catchment. This means that the improved model can increase the flood simulation accuracy on the slope by affording an accurate roughness on the slope.

The improved and raw model performance indices are shown in Table 3. The indices show that the improved model performs better than the raw model. The indices did not markedly increase in the large basin, but they obviously increased in a small catchment.

Figure 11 shows the spatial distributions of the flood peak anomaly percentage during 1998 to 1999 in the Yellow River basin. The FPAP is less than zero at the slope because the dynamic Manning's n is greater than the static Manning's n at the slope. However the FPAP is close to zero in the main river channel because the dynamic Manning's n is less than the static Manning's n in the river channel. The minimum FPAP is −45%, which shows that the simulated flood peak may have a large error at the slope with the static Manning's n, although the model has high performance at the basin outlet (Table 3).

The hourly flood peak is larger than the daily flood peak at the same time, and the hourly low flow is smaller than the daily low flow. The hourly n will be more variable than the daily n, resulting in a variable flow cross section (Equation (1)). Therefore, the hourly performance will be better than the daily performance using the dynamic n model. However, it is very difficult to obtain hourly discharge data over a long period. We will continue to study the hourly data in the future.

Figure 12 shows the seasonal spatial distribution of n in 1998 for the Yellow River. The n values were calculated by Equation (2) in the DTVGM. The n values are very small along the river channel and large on the slopes. The main causes are the higher vegetation density and shallower water depth on the slopes. The n values are small in spring and summer in the river channel because the water is deep during those two seasons. However, n is large in spring and summer on the slopes because the vegetation is luxuriant in spring and summer. The reverse result is shown in autumn and winter.

There are three advantages for applying the dynamic n in the hydrologic model. (1) The dynamic n is close to the true n value, and it is better than a constant in the raw hydrologic model. A more physical mechanistic parameter can decrease the model uncertainty and sensitivity of the parameter. (2) The improved simulated discharge will have a higher accuracy than the raw model in sub-basins. (3) The satisfactory spatial-temporal distribution of n can improve the flash flood and low flow simulations.

The weaknesses of the model are as follows: (1) the dynamic n equation is a statistic formula, making it difficult to calibrate the parameters; and (2) higher resolution DEM, vegetation and soil structure datasets are needed to calculate the dynamic n.

This research proposed a new way to provide a dynamic spatial-temporal distribution of Manning's roughness coefficients for hydrological models in basins. Based on the new concept, a distributed hydrological model was improved and applied to the Yellow River in China. The results have already shown that the improved model can provide more accurate simulated streamflow than the raw model. The new scheme of n took the LAI, soil components, hydraulic radius and water flow cross-section area into consideration, producing new n values that are closer to the true values than those from the raw model, especially on slopes. The n values are very small along the river channel and large on the slopes. The simulated flood peak may have a large error at the slope with the static Manning's n, although the model has high performance at the basin outlet. The simulated streamflow with the dynamic Manning's n is more accurate than that with the static Manning's n. The dynamic spatial-temporal distribution of n can now be used in other hydrological models or land surface models.

In the future, we will keep working to improve the roughness coefficient equation based on higher resolution data and other observed information. Additionally, we will develop a global routing in the land surface model using the new n scheme.

This study was supported by the Natural Science Foundation of China (No. 41475093), the Intergovernmental Key International S&T Innovation Cooperation Program (No. 2016YFE0102400) and the State Key Laboratory of Earth Surface Processes and Resource Ecology Open Research Program (No. 2017-KF-17).