Adaptive Determination of the Flow Accumulation Threshold for Extracting Drainage Networks from DEMs

: Selecting the ﬂow accumulation threshold (FAT) plays a central role in extracting drainage networks from Digital Elevation Models (DEMs). This work presents the MR-AP (Multiple Regression and Adaptive Power) method for choosing suitable FAT when extracting drainage from DEMs. This work employs 36 sample sub-basins in Hubei (China) province. Firstly, topography, the normalized difference vegetation index (NDVI), and water storage change are used in building multiple regression models to calculate the drainage length. Power functions are ﬁt to calculate the FAT of each sub-basin. Nine randomly chosen regions served as test sub-basins. The results show that: (1) water storage change and NDVI have high correlation with the drainage length, and the coefﬁcient of determination (R 2 ) ranges between 0.85 and 0.87; (2) the drainage length obtained from the Multiple Regression model using water storage change, NDVI, and topography as inﬂuence factors is similar to the actual drainage length, featuring a coefﬁcient of determination (R 2 ) equal to 0.714; (3) the MR-AP method calculates suitable FATs for each sub-basin in Hubei province, with a drainage length error equal to 5.13%. Moreover, drainage network extraction by the MR-AP method mainly depends on the water storage change and the NDVI, thus being consistent with the regional water-resources change.


Introduction
Water is a crucial resource for humans and the environment. Its distribution and flow characteristics have unique geographical traits. Extracting hydrologic information based on DEMs has received substantial attention in the hydrologic analysis of river networks. Classic problems in distributed hydrological modeling are the calculation of flow accumulation and the extraction of drainage networks [1,2]. A key variable in the extraction of drainage networks from DEMs is the flow accumulation threshold (FAT) [3]. The flow accumulation threshold is related to the upslope contributing area at a given point in a drainage area represented by a DEM. The FAT represents the minimum number of upslope cells flowing into a flow-receiving cell that is part of a drainage network (notice the FAT is a dimensionless variable). Applying an FAT to flow accumulation values delineates the drainage network from a DEM.
Several algorithms have been proposed for the extraction of drainage networks from DEMs, such as that of Callaghan and Mark, which applied the FAT as the basis for delineating drainage networks [4]. Several authors treat the FAT as a constant; Camara reported that the drainage networks obtained from DEM represent intermittent watercourses when the FAT equaled 20,000 m 2 (the areas of eight study sites equaled 1.64 km 2 ) [5]. This finding producing runoff that carves the river network through long-term erosion and deposition, while dense vegetation influences erosion [28,29]. This study's data stem from 36 regions employed for model development and 9 test areas in Hubei province, which are shown in Figure 1. These regions are sub-basins each with an area of more than 400 km 2 and they encompass most of the western and northern portions of Hubei province.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 19 by climate, vegetation, and topography. Mountains generate convection and frontal heavy rain producing runoff that carves the river network through long-term erosion and deposition, while dense vegetation influences erosion [28,29]. This study's data stem from 36 regions employed for model development and 9 test areas in Hubei province, which are shown in Figure 1. These regions are sub-basins each with an area of more than 400 km 2 and they encompass most of the western and northern portions of Hubei province.

Topographic Data
The digital elevation model (DEM) of the study area was collected from the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (Aster GDEM), having a 30-m spatial resolution. The Version 3 data have employed cloud masking to avoid cloud-contaminated pixels [30], which were downloaded from the Earthdata website (https://doi.org/10.5067/ASTER/ASTGTM.003, accessed on 21 November 2019).
The surface roughness (r) of the study region is defined as the ratio of the area of the actual surface to that of a smooth surface having the same geometric shape and dimensions. Equations (1) and (2) yield the surface roughness [31]: where SOD denotes the slope of the DEM's grid expressed as an angle in radians, geometric surface denotes the smooth surface area of the gird, and actual surface means the actual surface area of the grid.

Vegetation Index Data
The Normalized Difference Vegetation Index (NDVI) from June through August 2018 was determined from the Moderate-resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13A3) Version 6 data, which were downloaded from the Earthdata website (https://doi.org/10.5067/MODIS/MOD13A3.006, accessed on 21 November 2019). These data complements National Oceanic and Atmospheric Administration (NOAA)'s

Topographic Data
The digital elevation model (DEM) of the study area was collected from the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (Aster GDEM), having a 30-m spatial resolution. The Version 3 data have employed cloud masking to avoid cloud-contaminated pixels [30], which were downloaded from the Earthdata website (https://doi.org/10.5067/ASTER/ASTGTM.003, accessed on 21 November 2019).
The surface roughness (r) of the study region is defined as the ratio of the area of the actual surface to that of a smooth surface having the same geometric shape and dimensions. Equations (1) and (2) yield the surface roughness [31]: where SOD denotes the slope of the DEM's grid expressed as an angle in radians, geometric surface denotes the smooth surface area of the gird, and actual surface means the actual surface area of the grid.

Vegetation Index Data
The Normalized Difference Vegetation Index (NDVI) from June through August 2018 was determined from the Moderate-resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13A3) Version 6 data, which were downloaded from the Earthdata website (https://doi.org/10.5067/MODIS/MOD13A3.006, accessed on 21 November 2019). These data complements National Oceanic and Atmospheric Administration (NOAA)'s Advanced Very High Resolution Radiometer (AVHRR) NDVI products and provides continuity for a time series' historical applications [32]. The NDVI were resampled to have the same spatial resolution as the DEM with the IDW (Inverse Distance Weighted) method.

Water Storage Change from Precipitation, Runoff, and Evapotranspiration
The GPM_3IMERGM 06 precipitation data with a spatial resolution of 0.1 • during the period of June to August in 2018 were obtained from the Goddard Earth Sciences Data and Information Services Center (GES DISC) website (Huffman, G.J., E.F. Stocker, D.T. Bolvin, E.J. Nelkin, Jackson Tan. https://doi.org/10.5067/GPM/IMERG/3B-MONTH/06, accessed on 21 November 2019). The Global Precipitation Measurement (GPM) satellite launched on 28 February 2014, which is the successor of the Tropical Rainfall Measuring Mission (TRMM) and has made three critical improvements, that is, improving the spatial resolution of precipitation products from 0.25 • into 0.1 • , providing more coverage (covering 60 • N−60 • S), and adding sensitivity to liquid and solid precipitation [33,34].
NASA's GlobFigureal Land Data Assimilation System (GLDAS) generates optimal fields of land surface states and fluxes by blending satellite-and ground-based observational data products using advanced land surface modeling and data assimilation techniques [35]. GLDAS-2.1 started on 1 January 2000 with a spatial resolution of 0.25 • . GLDAS_NOAH025_M_2.1 products (including Evapotranspiration (ET) and Storm surface runoff (R)) covering this study area were used in this study. All of the data were downloaded from the GES DISC website (Beaudoing, H. and M. Rodell. https://doi.org/10.506 7/SXAVCZFAQLNO, accessed on 21 November 2019) [36].
Firstly, the precipitation, surface runoff, and evapotranspiration were resampled to 30 m spatial resolution with the IDW method. Secondly, the change in water storage (∆S) in any time interval can be calculated by the precipitation (P), surface runoff (R), and evapotranspiration (ET) occurring in an interval i according to the following function [37]:

Actual Drainage Networks Data
The OpenStreetMap (OSM) is a free map that can be created, edited, and updated by volunteers globally. In some urban areas, the quantity of networks already surpasses that of commercial or governmental datasets [38]. The actual drainage networks were obtained from the gis_osm_waterways_free_1 datasets corresponding to China. All of the data were downloaded from the OpenStreetMap website (https://download.geofabrik. de/asia/china-latest-free.shp.zip, accessed on 21 November 2019). They were employed in this work in the determination of the length of the actual drainage networks and the optimal determination of the FAT.

Drainage Networks Extraction
This paper's method for drainage network extraction applies the Deterministic eightnode (D8) algorithm of the hydrology module of ArcGIS, which consists of the following four steps: pit-filling, calculation of flow direction, calculation of flow accumulation grids, and drainage network extraction [4]. The specification of the FAT is required for drainage network extraction.
The D8 algorithm specifies each cell has its own flow direction, and the extraction of the drainage network conforms to this specification. This paper employs catchments whose drainage networks are linear and applies the drainage length method for the determination of the optimal FAT. The drainage length equation is calculated as follows: where L denotes the drainage length and l i represents the length of the i-th channel, where i= 1, 2, . . . , N [34].

Multiple Regression and Adaptive Power (MR-AP) Method
The Multiple Regression and Adaptive Power (MR-AP) method herein proposed resorts to a sub-basin's topography, the NDVI, and water storage change (∆S) as explanatory variables to accurately identify the FAT. The MR-AP method features two parts: multiple regression fitting to estimate drainage network length, and fitting of a power function between the FAT and drainage length to each sub-basin. The drainage network extracted with the MR-AP method is called the MR-AP network.
A workflow diagram of the MR-AP method is displayed in Figure 2. where L denotes the drainage length and li represents the length of the i-th channel, where i= 1, 2, …, N [34].

Multiple Regression and Adaptive Power (MR-AP) Method
The Multiple Regression and Adaptive Power (MR-AP) method herein proposed resorts to a sub-basin's topography, the NDVI, and water storage change (∆S) as explanatory variables to accurately identify the FAT. The MR-AP method features two parts: multiple regression fitting to estimate drainage network length, and fitting of a power function between the FAT and drainage length to each sub-basin. The drainage network extracted with the MR-AP method is called the MR-AP network.
A workflow diagram of the MR-AP method is displayed in Figure 2.

Multiple Stepwise Regression Fitting
The rationale for employing topography, NDVI, and water storage change as explanatory variables in drainage network extraction is founded on physical principles. Surface runoff derives from precipitation: the larger the precipitation, the larger the runoff, which has erosive power whether as overland flow or channelized [39]. Erosion is the beginning stage of drainage network formation. The topography is a key erosion factor: terrain with small surface roughness is less susceptible to erosion by runoff than terrain with large surface roughness, all other factors being equal. Vegetation, including forests and grasslands, on the other hand, increases the surficial shear capacity of soil or rock in which its roots are embedded and provides cover to the soil against a detachment of soil particles by falling rainfall. In other words, vegetation reduces the erosive power of runoff, besides increasing the water holding capacity of the soil and rock. These considerations motivate an exploratory analysis to ascertain the predictive skill of explanatory variables for the drainage length. This study applies the Pearson correlation analysis to test the correlation between drainage length and explanatory variables. Multiple regression employs topography, NDVI, and water storage change as the explanatory variables of drainage length of the 36 sample catchments. The fitted multiple regression is presented in the Results section.

Power Function Fitting for Each Sub-Basin
The FAT values selected as locally suitable decrease with increasing drainage length. Several authors have reported the FAT and drainage length are well related through power functions [14,40]. Likewise, this paper fits power functions in each sub-basin, as expressed by Equation (5): where L means the drainage length, a, b, and α denote the coefficients of the power function, and FAT denotes the flow accumulation threshold. A fraction equal to 1% of the basin area was applied as the initial FAT. The D8 algorithm was employed to calculate the drainage length corresponding to each initial FAT. Equation (3) of each sub-basin was fitted with the initial FAT and its corresponding drainage length. Next, the drainage length calculated by MR was applied to calculate the optimal FAT. This was followed by the extraction of the drainage network from the DEM in each sub-basin based on the optimal FAT. The drainage network extracted in this way herein is called the MR-AP drainage network.

Validation
The MR-AP method builds the multiple linear regression model based on the 36 sample sub-basins, then fits power functions to each sub-basin's data. The MR-AP method was subsequently evaluated with nine test sub-basins. The drainage length error was calculated by comparing extracted lengths with the actual drainage river length.
Analysts evaluate the quality of the extracted networks by comparing them with networks drawn on maps, represented in remotely sensed images, or represented in aerial photographs [23,41]. This paper compares the MR-AP drainage network with the actual drainage networks (graphed in Section 4.3) and the remote sensing image (graphed in Section 4.3) in evaluating the quality of the extracted network. The drainage length error is herein employed to assess the accuracy of the MR-AP drainage networks according to the following function:l wherel denotes the drainage length error, L i represents the length of all the extracted river network in a sub-basin, and L denotes the length of the actual drainage network in a sub-basin. A positive (negative) value ofl represents an overestimation (underestimation) of the drainage length.

Results of Drainage Length Using Multiple Regression (MR) Analysis
Many previous works have evaluated the spatial dependence of environmental parameters on the drainage length [42,43]. This study investigated the potential of using the explanatory variables cited above in building the multiple linear regression model.
The calculated Pearson correlation between actual drainage length and explanatory variables is listed in Table 1. The Pearson correlation values of the water storage change (∆S), NDVI, and maximum surface roughness (R max ) are equal to 0.85, 0.874, and 0.425, respectively, which indicates that all explanatory variables were significantly correlated with the actual drainage length at the 0.01 significance level (two-sided). Multiple linear regression model was applied next to account for the interaction between explanatory variables in the search for optimal regressions predicting the drainage network length.  Table 2 shows results for the chosen independent variables of the MR model, namely, the maximum surface roughness (R max ), the NDVI, and the water storage change (∆S). The constant coefficient, the coefficients associated with the water storage change, NDVI, and the maximum surface roughness were equal to −4270.479, 9015.26, 107.426, and −888.753, respectively. The R 2 (coefficient of determination) and adjusted R 2 (see Equations (A1) and (A2) in Appendix A) were 0.714 and 0.688, which means a relatively high degree of statistical association. The multiple regression equation is as follows: where L denotes the fitted drainage length (m), ∆S represents water storage change (mm), NDVI denotes the average NDVI, and R max is the maximum surface roughness of the sub-basins. The drainage length calculated in Equation (7) is called the MR drainage network length.

Drainage Network Results Using Adaptive Power (AP) Analysis
The drainage length and the FATs are well related by a power function law in which the power coefficient is negative [40]. For this reason, it has been argued the optimal threshold is defined by the point at which the second derivative of a power-law fitted to the drainage length as a function of the FAT approaches zero. This paper analyzed the 45 sub-basins (including 36 sample sub-basins and 9 test sub-basins) by fitting the FATs against the drainage with power functions. The drainage lengths of all the sub-basins are obtained with the MR method and the corresponding fitted FATs were calculated subsequently. Figure 3 shows the R 2 of all the sub-basins including the power functions corresponding to the best (No. 1) and worst (No. 36) sub-basins. All the sub-basins exhibit power functions between the FAT and drainage length: the range of the R 2 varied from 0.9341 to 0.9967. Figure 4 shows that the sub-basin No. 1 had the best R 2 , slightly lower than 0.997 (Figure 4a), and the sub-basin No. 36 has the lowest R 2 , slightly higher than 0.93 (Figure 4b). The average value and the standard deviation of R 2 of the power functions (see Equations (A3) and (A4) in Appendix A) equal 0.9795 and 0.0152, respectively, which means that the fitting of the function is accurate. It is herein established, therefore, that the FAT may be accurately estimated in terms of the calculated drainage length. sub-basins (including 36 sample sub-basins and 9 test sub-basins) by fitting the FATs against the drainage with power functions. The drainage lengths of all the sub-basins are obtained with the MR method and the corresponding fitted FATs were calculated subsequently. Figure 3 shows the R 2 of all the sub-basins including the power functions corresponding to the best (No. 1) and worst (No. 36) sub-basins. All the sub-basins exhibit power functions between the FAT and drainage length: the range of the R 2 varied from 0.9341 to 0.9967. Figure 4 shows that the sub-basin No. 1 had the best R 2 , slightly lower than 0.997 (Figure 4a), and the sub-basin No. 36 has the lowest R 2 , slightly higher than 0.93 ( Figure  4b). The average value and the standard deviation of R 2 of the power functions (see Equations (A3) and (A4) in Appendix A) equal 0.9795 and 0.0152, respectively, which means that the fitting of the function is accurate. It is herein established, therefore, that the FAT may be accurately estimated in terms of the calculated drainage length.  sub-basins (including 36 sample sub-basins and 9 test sub-basins) by fitting the FATs against the drainage with power functions. The drainage lengths of all the sub-basins are obtained with the MR method and the corresponding fitted FATs were calculated subsequently. Figure 3 shows the R 2 of all the sub-basins including the power functions corresponding to the best (No. 1) and worst (No. 36) sub-basins. All the sub-basins exhibit power functions between the FAT and drainage length: the range of the R 2 varied from 0.9341 to 0.9967. Figure 4 shows that the sub-basin No. 1 had the best R 2 , slightly lower than 0.997 (Figure 4a), and the sub-basin No. 36 has the lowest R 2 , slightly higher than 0.93 ( Figure  4b). The average value and the standard deviation of R 2 of the power functions (see Equations (A3) and (A4) in Appendix A) equal 0.9795 and 0.0152, respectively, which means that the fitting of the function is accurate. It is herein established, therefore, that the FAT may be accurately estimated in terms of the calculated drainage length.

Validation on the MR-AP Drainage Networks
To further verify the reliability of the MR-AP method in quantitative and qualitative ways, the actual drainage network data from OpenStreetMap were employed as the reference to quantitatively evaluate the quality of the MR-AP drainage network (see Figures 5 and 6), and to identify the exact location of the drainage network, as shown Remote Sens. 2021, 13, 2024 9 of 18 in Figure 7. Nine test sub-basins within Hubei province were randomly chosen for this evaluation, and they are depicted in Figure 1.

Validation on the MR-AP Drainage Networks
To further verify the reliability of the MR-AP method in quantitative and qualitative ways, the actual drainage network data from OpenStreetMap were employed as the reference to quantitatively evaluate the quality of the MR-AP drainage network (see Figures  5 and 6), and to identify the exact location of the drainage network, as shown in Figure 7. Nine test sub-basins within Hubei province were randomly chosen for this evaluation, and they are depicted in Figure 1. The MR-AP drainage network is calculated with the adaptive method depicted in Figure 2. Figure 5 compares the drainage length extracted by the MR-AP and actual drainage length, where La represents the actual drainage length and Le represents the MR-AP drainage length. The regression function is well fitted with R 2 (see Eqaution (A1) in Appendix A) equal to 0.9901 and correlation coefficient (CC, see Eqaution (A5) in Appendix A) equal to 0.995, which are slightly lower than 1. The root mean square error (RMSE, see Eqaution (A6) in Appendix A) equals 4736 m and the maximum absolute deviation (MAD, see Eqaution (A7) in Appendix A) equals 4250 m, while the minimum length of the drain-

Validation on the MR-AP Drainage Networks
To further verify the reliability of the MR-AP method in quantitative and qualitative ways, the actual drainage network data from OpenStreetMap were employed as the reference to quantitatively evaluate the quality of the MR-AP drainage network (see Figures  5 and 6), and to identify the exact location of the drainage network, as shown in Figure 7. Nine test sub-basins within Hubei province were randomly chosen for this evaluation, and they are depicted in Figure 1. The MR-AP drainage network is calculated with the adaptive method depicted in Figure 2. Figure 5 compares the drainage length extracted by the MR-AP and actual drainage length, where La represents the actual drainage length and Le represents the MR-AP drainage length. The regression function is well fitted with R 2 (see Eqaution (A1) in Appendix A) equal to 0.9901 and correlation coefficient (CC, see Eqaution (A5) in Appendix A) equal to 0.995, which are slightly lower than 1. The root mean square error (RMSE, see Eqaution (A6) in Appendix A) equals 4736 m and the maximum absolute deviation (MAD, see Eqaution (A7) in Appendix A) equals 4250 m, while the minimum length of the drain- The MR-AP drainage network is calculated with the adaptive method depicted in Figure 2. Figure 5 compares the drainage length extracted by the MR-AP and actual drainage length, where La represents the actual drainage length and Le represents the MR-AP drainage length. The regression function is well fitted with R 2 (see Eqaution (A1) in Appendix A) equal to 0.9901 and correlation coefficient (CC, see Eqaution (A5) in Appendix A) equal to 0.995, which are slightly lower than 1. The root mean square error (RMSE, see Eqaution (A6) in Appendix A) equals 4736 m and the maximum absolute deviation (MAD, see Eqaution (A7) in Appendix A) equals 4250 m, while the minimum length of the drainage network among the test basins equals 51,231 m. Figure 6 depicts the degree of normality of the regression standardized residuals. In terms of the histogram (Figure 6a), the mean of the regression standardized residual is −4.19 × 10 −15 , compared with 0 representing the best fitting, and the standard deviation equals 0.935, compared with 1 representing the best fitting. The P-P Plot (Figure 6b)    The similarity between the MR-AP extracted and the actual drainage networks for the nine test sub-basins is quantified with the absolute drainage length error listed in Table 3 error ranged between −6.25% and 9.87%, with the average drainage length error equal to 0.051. A negative drainage length error means an overestimation of the actual drainage by the MR-AP model, and vice versa. The miss hits, false hits, and true hits were calculated to evaluate the MR-AP drainage against actual drainage, and are listed in Table 3. The test sub-basins 3 and 8 exhibit the best and worst performances, respectively, as can be seen in Figures 7 and 8.    Figures 7 and 8 shows a qualitative comparison of the MR-AP drainage networks with the actual drainage networks and the remote sensing image. The actual drainage networks may have disconnected networks in some sub-basins, most of which cannot be captured by the MR-AP drainage network. The disconnected networks are always caused by human-made action, such as by canals and ditches in test 6 and test 8 in Figure 7. The OpenStreetMap drainage networks, which are created and updated by the field investigation of volunteers, may not be accurate in rough terrain and inaccessible areas. On the other hand, the MR-AP drainage networks extracted from DEMs that are generated by remote sensing methods are not limited by such terrain traits, as shown by test 3, 4, and 7 displayed in Figure 7. Despite these differences, the actual and the MR-AP drainage networks exhibit good resemblance, such as demonstrated in test basin 3 displayed in Figure 7. The MR-AP drainage networks are compared with remote sensing imagery from Google Earth Pro in Figure 8 [44]. Gullies, which are generated by runoff erosion, promote the convergence of runoff. The distribution of gullies and surface runoff exhibited a good degree of consistency, and for this reason, gullies are considered incipient drainage networks. Gullies are well represented in the MR-AP drainage networks depicted in Figure 8. More detailed inspection of the extracted river networks in test sub-basins 3, 4, and 7 can be garnered from Figure 9. These results establish the MR-AP method is a suitable estimator of the FAT and an accurate tool for extracting accurate drainage networks. These results establish the MR-AP method is a suitable estimator of the FAT and an accurate tool for extracting accurate drainage networks.

The Statistical Association between the Drainage Length and the Explanatory Variables
The test results presented in Section 4.3 establish that the calculated results based on explanatory variables (water storage change, NDVI, and maximum surface roughness) performed well compared with the OSM drainage networks, which demonstrated that those explanatory variables can be effective environmental variables for extracting drainage networks with the MR-AP method. The statistical association between the drainage length and the explanatory variables in the 45 areas was further explored. Figure 10 shows the statistical association between the drainage length derived from the OSM drainage network and the explanatory variables. It can be seen in Figure 10 that there is a significant positive statistical association between the drainage length and the water storage change and the NDVI, and there is a high correlation in these regressions. However, it is also evident in Figure 10 that there is a weak statistical association between the maximum surface roughness and the drainage length; yet, several studies have reported that the surface roughness has an effect on the drainage length [15,21]. The weak statistical association between the maximum surface roughness and the drainage length displayed in Figure 10 may be due to the fact that some of the test sub-basin (such as test 6 and test 8) are located in plain areas, and the DEM-based drainage network extraction method does not perform well in flat terrain [45,46].

Comparison between MR-AP, Mean Change-Point Analysis, and the Fitness Index Method
Change-point analysis is applied to detect an abrupt or structural change in the distributional properties of data [47], and it has been used to calculate the optimal threshold in discrete models in many studies [48]. The fitness index method is also used to calculate the optimal threshold in drainage network extraction [49]. Figure 11 depicts the drainage length error, true hits, miss hits, and false hits (see Equations (A8), (A9) and (A10) in Appendix A) of the extracted drainage networks obtained with the MR-AP model, mean change-point analysis, and the fitness index method. It can be seen in Figure 11a that the drainage extracted by MR-AP consistently features the best performance with respect to

Comparison between MR-AP, Mean Change-Point Analysis, and the Fitness Index Method
Change-point analysis is applied to detect an abrupt or structural change in the distributional properties of data [47], and it has been used to calculate the optimal threshold in discrete models in many studies [48]. The fitness index method is also used to calculate the optimal threshold in drainage network extraction [49]. Figure 11 depicts the drainage length error, true hits, miss hits, and false hits (see Equations (A8)-(A10) in Appendix A) of the extracted drainage networks obtained with the MR-AP model, mean change-point analysis, and the fitness index method. It can be seen in Figure 11a that the drainage extracted by MR-AP consistently features the best performance with respect to the drainage length error, with the absolute error under 10%. The fitness index method tends to underestimate the drainage length, having negative errors. The three methods exhibit opposite quality of performance with respect to the miss hits and false hits according to Figure 11b,c. For example, the MR-AP model exhibits the best performance with respect to miss hits and the worst performance with respect to false hits in many sub-basins. Figure 11d displays the true hits obtained with the three methods; the MR-AP method features the best true hits in eight sub-basins and worse performance than the mean change-point analysis only in test sub-basin 6. These results establish that the proposed MR-AP model extracts drainage networks with higher accuracy than the other two methods.

Conclusions
This work chose 45 sub-basins in Hubei province and evaluated three potential explanatory variables for extracting drainage networks, namely, maximum surface roughness, NDVI, and water storage change. This study proposes the MR-AP method to calculate the FAT adaptively, which employs power functions relating the FAT and drainage length in each sub-basin. This paper's results indicate that (1) it is reasonable and accurate to employ maximum surface roughness, normalized difference vegetation index, and water storage change as explanatory variables for the drainage network length; (2) multiple regression fitting provides superior fitting of drainage network length compared to individual regressions between the drainage length and each explanatory variable; (3) the drainage length and the FATs exhibited strong correlation expressed through power functions fitted with data for Hubei province, China; (4) the MR-AP method was shown to be accurate in estimating the FAT and for extracting drainage networks. However, drainage networks are created by complex processes that involve a number of factors. Future re-

Conclusions
This work chose 45 sub-basins in Hubei province and evaluated three potential explanatory variables for extracting drainage networks, namely, maximum surface roughness, NDVI, and water storage change. This study proposes the MR-AP method to calculate the FAT adaptively, which employs power functions relating the FAT and drainage length in each sub-basin. This paper's results indicate that (1) it is reasonable and accurate to employ maximum surface roughness, normalized difference vegetation index, and water storage change as explanatory variables for the drainage network length; (2) multiple regression fitting provides superior fitting of drainage network length compared to individual regressions between the drainage length and each explanatory variable; (3) the drainage length and the FATs exhibited strong correlation expressed through power functions fitted with data for Hubei province, China; (4) the MR-AP method was shown to be accurate in estimating the FAT and for extracting drainage networks. However, drainage networks are created by complex processes that involve a number of factors. Future research should consider other factors such as temperature, humidity, soil type, and surface geologic features as possible explanatory variables in the extraction of drainage networks. Multiple regression models require a linear statistical association between the dependent and explanatory variables. Nonlinear fitting models could also be considered in future studies. The morphology of extracted drainage networks depends on the calculation method for flow direction. Multi-flow algorithms such as the D-infinity method could be considered instead of the single-flow D8 algorithm in future works.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In Table 2, the R 2 and Modified R 2 represents the coefficient of determination and the adjusted coefficient of determination of the MR model, respectively, which are calculated with the following equations: where L, L, andL denote the actual, average, and predicted drainage length, respectively, N denotes the number of the sub-basins, and k denotes the number of independent variables. In Figure 3 the Mean and StdDev represents the mean and the standard deviation of the R 2 , respectively, which are calculated with the following equations: where R 2 i denotes the R 2 of the i-th sub-basin and N denotes the total number of the sub-basins.
In Figure 5 La represents the actual drainage length (La, m) and Le represents the drainage length extracted by MR-AP (Le, m). CC, RMSE, and MAD represent the correlation coefficient, root mean squared error, and median absolute deviation between La and Le, respectively, which are calculated with the following equations: where La i and Le i denote the actual drainage length and the drainage length extracted with the MR-AP method in the i-th test sub-basin, respectively, La and Le denote the average La and Le in all test sub-basins, and n denotes the total number of test sub-basins.
In Table 3 and Figure 11, the true hits, miss hits, and false hits represent the truly or correctly extracted, mis-extracted or failed to be extracted, and falsely or incorrectly extracted drainage networks when compared to the actual drainage networks, respectively, which are calculated with the following equations: where L, L true , andL denote the actual, true, and predicted drainage lengths respectively.