Basin-scale high-resolution extraction of drainage networks using 10-m Sentinel-2 imagery

Extraction of drainage networks is an important element of river flow routing in hydrology and large-scale estimates of river behaviors in Earth sciences. Emerging studies with a focus on greenhouse gases reveal that small rivers can contribute to more than half of the global carbon emissions from inland waters (including lakes and wetlands). However, large-scale extraction of drainage networks is constrained by the coarse resolution of observational data and models, which hinders assessments of terrestrial hydrological and biogeochemical cycles. Recognizing that Sentinel-2 satellite can detect surface water up to a 10-m resolution over large scales, we propose a new method named Remote Sensing Stream Burning (RSSB) to integrate high-resolution observational flow location with coarse topography to improve the extraction of drainage network. In RSSB, satellite-derived input is integrated in a spatially continuous manner, producing a quasi-bathymetry map where relative relief is enforced, enabling a fine-grained, accurate, and multitemporal extraction of drainage network. RSSB was applied to the Lancang-Mekong River basin to derive a 10-m resolution drainage network, with a significant reduction in location errors as validated by the river centerline measurements. The high-resolution extraction resulted in a realistic representation of meanders and detailed network connections. Further, RSSB enabled a multitemporal extraction of river networks during wet/dry seasons and before/after the formation of new channels. The proposed method is fully automated, meaning that the network extraction preserves basin-wide connectivity without requiring any postprocessing, hence facilitating the construction of drainage networks data with openly accessible imagery. The RSSB method provides a basis for the accurate representation of drainage networks that maintains channel connectivity, allows a more realistic inclusion of small rivers and streams, and enables a greater understanding of complex but active exchange between inland water and other related Earth system components.


Introduction
Extraction of drainage networks facilitates answering fundamental questions about how water flows on the Earth's surface.The extraction of drainage networks has been traditionally rooted in hydrologic science where it is particularly applied in river flow routing (David et al. 2011;P. Lin et al., 2018b;Yamazaki et al. 2013).However, high-resolution drainage networks have been increasingly used to support large-scale hydrologic predictions and spatially detailed research applications including the assessments of flood inundation, dam failure, and reservoir operation (Lehner and Grill 2013;Shin et al. 2020;Yamazaki et al. 2019).
The role of smaller streams within river networks in global biogeochemical cycling remains poorly understood and is likely underestimated (Benstead and Leigh 2012;Raymond et al. 2013).Traditional approaches assume that rivers and streams only transport materials from the land to the ocean.Recent evidence, however, has gradually unveiled the considerable amounts of materials exchanged and stored in river networks (Aufdenkampe et al. 2011;Bastviken et al. 2011;Battin et al. 2009Battin et al. , 2008;;Butman and Raymond 2011;Cole et al. 2007), indicating that the river system is an essential component of global biogeochemical cycles that have been long overlooked.Within river networks, small rivers are more biogeochemically active, partly due to their high density along with their intense interactions with the atmosphere and benthic substrate, which contribute to their significant total impact (Butman et al. 2016;Butman and Raymond 2011;Raymond et al. 2013).For instance, small rivers (i.e., those corresponding to Horton-Strahler stream order (Strahler 1957) of 3 or less) account for ~89% of the global river length, emitting more than half of the greenhouse gases from the Earth's freshwater including lakes and wetlands (Butman et al. 2016;Raymond et al. 2013), even though they represent a relatively small portion of surface area (Downing et al. 2012).Small rivers exhibit high evasion rates of greenhouse gases, often two to three times those of larger rivers (Aufdenkampe et al. 2011), implying that small discrepancies in the estimates of small rivers can incur large uncertainties in global estimates of greenhouse gas emissions.Another high uncertainty arises from the poorly known temporal patterns of small rivers (Costigan et al. 2015;Stanley et al. 1997).The expansion and contraction of rivers over time can result in significant biogeochemical effects (Butman et al. 2016;HILL et al. 2010;Hotchkiss et al. 2015;Marx et al. 2017), and such changes are relatively extreme for small rivers (Allen et al. 2018).
Delineation of smaller streams and rivers depends upon a higher resolution of models or maps because rivers and streams are often delineated at large scales by either models or maps.With the advent of the digital elevation model (DEM), automatic extraction of drainage networks has become convenient (Tarboton et al. 1991) with the pixel-wise, directionbased methods, such as D8 (Ocallaghan and Mark 1984) and D-infinity (Tarboton 1997).The recent emergence of high-resolution DEM (<10-m), mostly captured by the light detection and ranging (LiDAR) technique, has increased the resolution of extracted rivers (Clubb et al. 2014;Hooshyar et al. 2016;Passalacqua et al. 2012).However, increased resolution also introduces complications in topography (e.g., low-relief landscape, vegetation, and artificial constructions), requiring substantially more effort and computations to distinguish river channels and construct river networks (Passalacqua et al. 2010;Yamazaki et al. 2019).Meanwhile, the additional processing of high-resolution DEM is often unable to preserve drainage networks across the entire landscape due to the focus on river channels only (Bryndal and Kroczak 2019;Wu et al. 2019), which largely limits hydrological routing analysis and scaling estimates of river networks for carbon budget.Most importantly, the availability of high-resolution DEM is now far from sufficient to cover the entirety of a large river basin (>10 5 km 2 ), preventing the extraction of drainage networks at large scales.
The map-based delineation of rivers typically relies on optical satellite imagery for major river basins or a large spatial extent.Although the Synthetic Aperture Radar (SAR) is an effective microwave-based imaging tool for flood inundation detection (Hostache et al. 2018), the complexity of processing SAR backscatter signals impedes its use in large-scale water mapping that still relies on optical imagery (Cian et al. 2018;Huang et al. 2018).Water extents are photogrammetrically mapped from remote sensing imagery by either supervised or unsupervised methods (Ozesmi and Bauer 2002).For an entire basin, unsupervised water index methods based on optical imagery are widely adopted owing to their high efficiencies (Huang et al. 2018), such as the Normalized Difference Water Index (NDWI) (McFeeters 1996), Modified Normalized Difference Water Index (MNDWI) (Xu 2006) and Multispectral Water Index (MuWI) (Wang et al. 2018).Based on water maps, river pixels and river centerlines can be further identified and measured (Pavelsky and Smith 2008).For example, Allen and Pavelsky (2018) updated the global extent of rivers and streams by applying the MNDWI method on Landsat images captured in a month of average river discharge, representing a state-of-the-art feat in the field of large-scale Z. Wang et al. river remote sensing.Moreover, multitemporal remote sensing images corresponding to multitemporal surface water extents can potentially be used to detect river dynamics (Barefoot et al. 2019;Donchyts et al. 2016a;Pekel et al. 2016).However, the major drawback of satellite imagery-based extraction is that it extracts many river centerlines rather than one river network per basin (Fig. 1).Although most extracted river centerlines are connected, the overall connectivity within a basin cannot be guaranteed and can easily be undermined by a single break point.Such a break point is often observed at a pixel corresponding to on-river objects (e.g., dams, bridges, ships) or water detection omissions (e.g., contamination by clouds, shadow, sunglint).This issue becomes more prominent as the resolution of the imagery increases and the spatial extent is enlarged.Lack of connectivity hinders the construction of the river networks.Even if the basin-wide connectivity can be theoretically reconstructed with some image processing techniques (e.g., Lin et al. 2010), no algorithm is yet available to construct the river networks solely from image-derived river centerlines due to difficulties in flow direction identification, not to mention to construct drainage networks for the entire globe.A second issue is the spatial resolution.It is challenging to directly detect rivers narrower than one-pixel size, which often corresponds to at least tens of meters for freely available imagery with large-scale geographical coverage, as more sophisticated subpixel analysis is uncommon in large-scale river networks detection.Meanwhile, the spatial resolution of the imagery is not high enough to extract small rivers partially covered by vegetation.Taller riparian vegetation with sufficient water availability worsen the coverage over the river (Dosskey et al. 2010;Mokgoebo et al. 2018).Last, a minor concern is that current surface water remote sensing at large scales is mostly based on Landsat imagery, which limits the river detection to a 30-m resolution.While the Landsat-like Sentinel-2 imagery has the potential to provide 10-m resolution for river delineation, which is particularly meaningful for small rivers as it can extend the river detection to smaller stream orders, Sentinel-2 imagery is still used far less in river studies compared to that of Landsat.
In contrast to the highly appreciated functions of small rivers within a river network, extraction of the detailed river networks over sufficiently large extents remains challenging (Benstead and Leigh 2012;Wu et al. 2019).Model-based extraction is necessary to infer small stream orders from the drainage networks of the entire landscape but is constrained by the lack of high-resolution data for large extents.Map-based extraction can provide observational river locations but lacks the basinwide connectivity necessary to construct the networks and is also limited by scarcity of higher resolution data at large scales.However, a comparison between the two types of data suggests that optical imagery is normally commercially cheaper than DEM at equivalent resolutions due to less intensive elaboration of the raw data (Grosse et al. 2012).It is also true for the openly accessible situation where Sentinel-2 can provide 10m resolution globally while global DEMs (e.g., AW3D, SRTM, ASTER, TanDEM-X) often have a spatial resolution of 30-m at best.
From this context, this study aims to investigate the integration of surface water remote sensing into drainage extraction over the entirety of a large river basin for a more realistic representation of small rivers.The primary goal is to examine whether optical remote sensing can be coupled with DEMs to improve the resolution and dynamic extraction of drainage networks, and, if so, explore how it can be integrated with DEMs.In that regard, we propose a new method rooted in the stream burning approach (Saunders 1999) with the input of a Sentinel-2-based water index MuWI (Wang et al., 2018), hereafter referred to as the Remote Sensing Stream Burning (RSSB).The Lancang-Mekong River basin was selected for testing the RSSB method because of its sufficiently large area and vastly diverse landscapes.To evaluate the spatial accuracy of our integrated method, we compared the extraction results with river centerlines (George H. Allen and Pavelsky, 2018).We then demonstrate the extraction of river dynamics by presenting the details of extracted river networks for multiple periods and over a meander.We have made the following advances over previous basin-to-global scale studies on drainage networks extraction: (1) higher resolution of drainage networks can be achieved at large scales based on the fusion of Sentinel-2 imagery and coarse-resolution DEM; (2) the increased resolution provides improvements in drainage networks accuracy with realistic meander representation and networks connections; (3) applying the proposed method to individual images enables the extraction of multitemporal drainage networks over different periods of interests, e.g., during wet/dry seasons and before/after the formation of new channels; (4) as a fully automated extraction is conducted for a large river basin, postprocessing that connects false endorheic watersheds is no longer required.Although we only implemented the framework using Sentinel-2 imagery over the Lancang-Mekong River basin, the new method presented here can potentially be incorporated into other high-resolution imagery and applied to other large river basins or even applied globally.This approach to extracting high-resolution drainage networks would provide a methodological basis for understanding the role of river networks, particularly that of smaller streams and rivers, in Earth system.

Study area
The Lancang-Mekong River basin (Fig. 1) is one of the major river basins in the world with transboundary coverage.The river basin is characterized by diverse fluvial geomorphology with valley-constrained regions upstream to the bedrock-constrained areas downstream (Meshkova and Carling 2012).The Lancang-Mekong River is the longest river in Southeast Asia, at ~5000 km in length (Liu et al. 2007), and has an reported annual discharge of ~145,000 m 3 /s (MRC 2010; Pokhrel et al. 2018).It covers an area greater than 795,000 km 2 and is home to more than 60 million people from six countries: China, Myanmar, Laos, Thailand, Vietnam, and Cambodia.The gradient of the upper Lancang River is approximately 2 m/km, more than ten times that of the lower Mekong River, indicating more convergent topography exists upstream while divergent but well-defined banks are prevalent downstream (Pokhrel et al. 2018).In the densely populated lower basin, the monsoon climate determines the high variability of surface water (Ruiz-Barradas and Nigam 2018; Yang et al. 2019).The boundary of the study area follows the basin boundary from the HydroBASINS dataset (Lehner and Grill 2013).

Reflectance data
The reflectance data used in this study were from a Sentinel-2 Multispectral Instrument (MSI) by the European Space Agency (ESA).The Sentinel-2 MSI provides 13 spectral bands at three spatial resolutions (10-m, 20-m, and 60-m) for the period of June 2015 and beyond.Before the Sentinel-2, Landsat images captured starting from January 1987 were used, from sensors of the Thematic Mapper (TM) on Landsat 5 and Operational Land Imager (OLI) on Landsat 8, all at the spatial resolution of 30 m.

Topography data
We used the Multi-Error-Removed Improved-Terrain DEM (MERIT DEM) at a resolution of 90-m (Yamazaki et al. 2017) as the topography data.We selected not to use the 30-m open access data (e.g., SRTM30) for two reasons: 1) we aimed to demonstrate the robustness of this method by using relatively lower-resolution topography; 2) MERIT DEM had multiple elevation errors (e.g., in mosaicking tiles) removed so that technical disturbance to the extraction of the drainage networks could be minimal (Supplementary Information).For the extent of the Lancang-Mekong River basin, MERIT DEM was processed from SRTM3 DEM that was originally acquired by an 11-day InSAR mission in February of 2000 (Farr et al. 2007).Multiple errors were corrected in the topography, including Z. Wang et al. absolute bias, stripe noise, speckle noise and tree height bias (Yamazaki et al. 2017), making it suitable for drainage networks extraction (Shin et al. 2020(Shin et al. , 2019) ) and for baseline comparison.Potential flow disconnectivity in extraction could be significantly mitigated from the edge effect caused by mosaicking inconsistent tiles (Lehner et al. 2008;Wu et al. 2008).It is worth noting that for water bodies, SRTM3 DEM (as well as the derived data) was intended to depict the elevation of the water surface at the time of the SRTM mission in February 2000 (Slater et al. 2006) instead of the bathymetry at other times.

Reference river location data
Reference river locations were obtained from the Global River widths from Landsat (GRWL) database (Allen & Pavelsky, 2018).GRWL was compiled from Landsat images in the near mean-discharge month circa 2015.The modified normalized difference water index (MNDWI) was applied in developing GRWL to identify water pixels from Landsat images.Further processing on extracted water mask eventually produced GRWL with information regarding the location, width, and braiding index of rivers, in both vector and raster forms.The river vector (polyline) and Z. Wang et al. river raster (one Landsat pixel wide) from GRWL were measured as the river centerline using RivWidth (Pavelsky and Smith 2008) software.Constrained by Landsat inputs, only rivers ≥30-m wide were included in GRWL, and the river connectivity was not guaranteed.We extracted river vectors from GRWL within our study area to obtain 649,189 river reaches; the average length of reach segmentation was 36.15 m

Methods
Three main steps were adopted for accurate and dynamic extraction of river networks (Fig. 2).First, remote sensing input was prepared by identifying water presence based on Sentinel-2 bands, and then MuWI percentiles were computed to create a composite image.Second, the water index was integrated with the coarse topography using the proposed Remote Sensing Stream Burning (RSSB) to inherit the highresolution from the imagery.Third, flow direction and flow accumulation were computed through the common hydrological routing for the connected drainage network.

Creating the composite image from reflectance percentile
We created the Sentinel-2 or Landsat composite by calculating perband, per-pixel percentiles of reflectance within a given period on the Google Earth Engine platform (Gorelick et al. 2017).An adequate percentile of reflectance over sufficient scenes can exclude the outlier values from clouds, shadows, and sunglint (Donchyts et al. 2016b(Donchyts et al. , 2016a;;Ortiz et al. 2017;Wang et al. 2018).We used the 50% percentile in this study because it can represent the average condition of a period while reducing uncertainties from intermittent acquisition time, such as variable water extent and vegetation disturbance (Donchyts et al. 2016b;Lin et al., 2018a;Ortiz et al. 2017).The exact percentile can be determined a priori using other datasets (e.g., Wilson and Jetz 2016).To produce the latest, most accurate Lancang-Mekong River network, we used a total of 13,376 Sentinel-2 scenes from 2015 to 06-23 to 2019-04-01 over the entire basin.High revisit frequency and overlap swath of the twin constellation ensures each pixel is contained within 50 to 157 scenes for percentile calculation.For the period before 2015, we used Landsat scenes over a selected region to demonstrate the dynamics of flow pathways.Additionally, to best eliminate cloud contamination, we filtered out low-quality pixels (e.g., pixels of opaque clouds, cirrus clouds, and cloud shadows) using quality-assessment bands (i.e., QA60 band of Sentinel-2 and pixel_qa band of Landsat).

Identifying water presence by multispectral water index (MuWI)
Water masks or water maps are usually produced using threshold values of water index with binary coding for water and non-water pixels (Allen and Pavelsky, 2018;Yamazaki et al. 2017).In this study, we extended the utility of water index by assuming its value was a relative quantification of surface water presence.This assumption could be reasonable because a lower percentage of water coverage corresponds to a lower value of water index at the pixel or subpixel level (Feyisa et al. 2014;Wang et al. 2018).
Among many water indexes, we selected the Multispectral Water Index (MuWI) to identify water presence because it is a native 10-m water index on Sentinel-2 with improved accuracy.Sentinel-2 has multiple band resolutions (10, 20, 60 m), indicating previous indexes deliver either lower resolution (e.g., MNDWI) or lower accuracy (e.g., NDWI).MuWI adapts to both Sentinel-2 and Landsat with improved low-albedo performance (Wang et al. 2018).Low-albedo pixels, originating from the shadows of buildings, hillslopes, and clouds, usually contribute to commission error and are speculated to be the primary error source in water mappings (Fisher et al. 2016).A systematic comparison shows that the commission error in water mappings using MuWI is ~6%, which is significantly lower than that from NDWI (~18%) and MNDWI (~14%) (Wang et al. 2018).We calculated MuWI based on the composite images (Eqs. 1 and 2 for Sentinel-2 and Landsat, respectively).Depending on the adopted satellite, MuWI provides 10-m water presence from Sentinel-2 and 30-m water presence from Landsat.
where ND(i,j) denotes the normalized difference of band i and band j; and ρ represents top-of-atmosphere (TOA) reflectance for corresponding Sentinel-2 or Landsat bands: (3)

Stream burning and flow routing
Stream burning is a widely adopted flow enforcement method to produce conditioned DEM for deriving drainage networks (Lehner and Grill, 2013;Saunders and Maidment 1996;Yamazaki et al., 2019).Conventional stream burning often uses a rasterized stream map as the mask to lower the stream pixel.
Taking this one step further, we incorporated MuWI-represented water presence as the burnt layer into the DEM to force the flow.A resolution of 10 m can be achieved when producing water presence using Sentinel-2 images.The representative moment, or period of water presence, is more flexible in selecting satellite images because the intermediate step of producing water maps (particularly, thresholding) is avoided.As MuWI-derived water presence layer has spatially continuous values for both water and land pixels, frequently reported issues resulting from only burning water pixels, such as disconnection (Turcotte et al. 2001;Yamazaki et al. 2019) and parallel channels (Yamazaki et al. 2015(Yamazaki et al. , 2012)), are expected to be mitigated.
where φ is the burning intensity parameter.Based on the conditioned DEM from stream burning, hydraulic conditioning was conducted to overcome negative effects induced from topographic depressions (Woodrow et al. 2016).We calculated flow direction and flow accumulation using the D8 algorithm (Ocallaghan and Mark 1984) primarily due to computational efficiency and consistency with previous datasets, especially in terms of large-scale applications.Following Lin et al. (2019), we applied a channelization threshold of 25 km 2 .Flow calculations were performed through the TauDEM software package (Tarboton 2005) under the Message Passing Interface (MPI) parallel computing mode on a 48-core workstation.

Analyses
We assessed the positional accuracy of the extracted river using Goodchild's measure for linear features (Goodchild and Hunter 1997).The proportion of extracted river that lies within the buffer of the reference river centerline was computed according to buffer sizes (Eq. 5 where L represents the length of the river vector line).With increasing buffer size, x, the positional difference between the extraction and the reference was estimated using the threshold p(x)≥0.85,as done in Donchyts et al. (2016b).

p(x) =
L(overlapped with the buffer of size x) L total (5) The statistical difference between the extraction and the reference was assessed using the two-sample Kolmogorov-Smirnov test (Allen et al. 2018;Venables and Ripley 2013).
Burning intensity parameter, φ, is the major controlling parameter for both topographical modification and the final extraction.We compared and contrasted various φ values from 1, 5, 10 and 20 m to characterize the sensitivity of the method.However, since the least modification was expected on the baseline DEM (Callow et al. 2007), we further analyzed the Pareto optimum where necessary burning intensity ensured replication of the river location, while further increasing Z. Wang et al. burning intensity only marginally improved accuracy.To quantify heterogeneous burning intensity, flow accumulation (A) either from D8 or D-infinity method is used as a basis: where A 95% represents the 95% percentile of flow accumulation in the entire basin.This equation can warrant smaller burning intensity for headwaters while larger burning intensity for downstream floodplains.

Higher-resolution extraction of the drainage network
To examine the possibility of increasing resolution of the drainage network, we used the RSSB method to extract the high-resolution drainage networks for the entirety of the selected sub-basin.Significant improvements are evident in the high-resolution extraction in terms of realistic river representations (Fig. 3).For comparison, equivalent Sentinel-2 imagery and DEM are integrated in both extractions, but the lower one adopts a high resolution (10-m) following the RSSB method while the upper one adopts a coarse-resolution (90-m) following the ordinary stream burning method.Not only does increased river-reach sinuosity indicate a higher level of fidelity, but an increase in authentic tributaries also suggests that improvements of networking correctness can be identified in the RSSB high-resolution extraction (Fig. 3a).Large structural errors in the coarse resolution (90-m) extraction are evident, likely due to low-relief terrain in the area and false tributary joints.Furthermore, the high-resolution extraction successfully excludes the oxbow lakes that are unilaterally connected to the river, which have completely different flow regimes, as shown in the yellow circled areas in Fig. 3a.At the networks region in Fig. 3a, an additional river length of 41% is delineated in the high-resolution extraction (57,460 m vs. 40,687 m), which is associated with a minimum of nearly half the difference in estimates of carbon emissions.A closer inspection at the meander level is shown in the example of Fig. 3b, where the meander is accurately delineated in the high-resolution extraction, but not in the ordinary coarse-resolution extraction.At this meander, ~100% longer river length is observed, implying that small rivers with more curvatures require such high-resolution extraction for more realistic estimates, particularly for biogeochemical processes.
The cumulative density functions (CDFs) statistically show that accuracy was improved in the RSSB-derived higher-resolution extraction (Fig. 4).The 90-m DEM-based 10-m extraction (red line) show the Table 1 shows the river length and river networks density distribution of drainage networks extracted from the two resolutions.For both extractions, the river length and networks density all roughly follow the classic power law (Leopold and Maddock 1953).However, the RSSBderived high-resolution drainage networks presents a longer river length, and thus, a denser network.This discrepancy can be in the order of 2.5 times for both total length and networks density.Notably, the high-resolution approach extracts one more stream order (ω = 9) compared to that of the coarse-resolution, implying that increased resolution not only improves representation of flowlines but also strengthens the level of details in the network.

Multitemporal representation of drainage network
The ability of RSSB to extract the dynamics of the river networks is vital.To demonstrate the multitemporal aspect of this method, we depict a reach in the Lancang-Mekong River basin where the major flow route is likely to change over time (Fig. 5a).In Area A, a new river channel initialed circa 2005, as can be inferred from annual MuWI change and supported by two images captured in 2001 and 2014 (Fig. 5b).In Area B, the water presence percentage varied significantly between dry and wet seasons in a U-turning channel area of the Lancang-Mekong River, i.e., February to May (dry) and June to October (wet) (MRC 2010), as shown in Fig. 5c.Four periods of wet/dry seasons before/after new-channel initiation differ from the water presence patterns as measured by MuWI (right panel of Fig. 5a).The variations of water presence propagate to the drainage networks extraction using the RSSB method.During the dry season after initiation of the new channel (green line in Fig. 5a), a higher water presence is apparent in the upper bifurcation near Area A, so that the flowline in the extraction of this period diverts into the new channel.Similarly, more water is present in the lower bifurcation near Area B during the wet season, after initiation of the new channel (red line in Fig. 5a).As a result, the extracted drainage flows through the lower meanders rather than upward through the U-shaped channel.All four extractions agree better with the river channels, whether bankfull or exposed, in contrast to the extraction performed without burning (gray line in Fig. 5a).Although the upper route is also connected during the wet season after the formation of the new channel, the single flow direction approach (D8) in our method identifies the lower route to be the major route.Despite possible debate on the criterion to determine the major route in this situation of multiple de facto flow directions, the results show that varying remote sensing inputs in the RSSB method could differentiate the flow routes during environmental changes (Fig. 5).As this multitemporal capability of RSSB method relies on the selection of imagery, it needs to be further tested and developed in future work (e.g., possible issues with mismatching scenes from a rapid avulsion).

Basin-wide extraction of the drainage network
To demonstrate the added spatial accuracy afforded by integrating surface water remote sensing, we extracted the drainage networks for the entire Lancang-Mekong River basin and compared it with the observational river centerlines.We note that due to computational constraints (>1 TB memory requirement), the basin-wide extraction is upscaled to 90 m from 10 m.Two extraction results (Origin and RSSB) were generated both at 90 m resolution; thus, the primary difference between them was whether or not Sentinel-2 imagery was integrated through MuWI.Even though it is not the best extraction promised by RSSB method, there is still a significant reduction in positional error after applying our method (Fig. 6).Overall, the distance to authentic river centerlines is reduced by 48%, from 236 m of the original

Table 1
Stream order ω, number n ω , mean length L mean , total length L total , and the river networks density D, of drainage networks from the high-resolution extraction using RSSB and the coarse-resolution extraction using ordinary stream burning.

High-resolution (RSSB)
Coarse-resolution (ordinary)   extraction to 123 m of the RSSB extraction.Wider reaches are associated with lower positional differences for both extraction results where river width is less than or equal to 500 m.However, the error-width relationship is not consistent for reaches wider than 500 m, likely because more bifurcations and braided sections exist in wider channels.All width categories show a 50±10% error reduction, suggestive of a systematic improvement from remote sensing treatment in terms of location accuracy.Spatial distribution of positional errors over the Lancang-Mekong River basin from the two extraction results is compared in Fig. 7.For both results, the positional error is higher in the lower basin than that in the upper basin (Figs.7a-b).This likely occurred because the flat and populated landscape, which hinders accurate extraction of the river network, is often observed in downstream regions of the Lancang-Mekong basin (Deng 2019;Viero et al. 2019).Fig. 7c reveals the improvement, or error reduction, of positional accuracy across most basin areas.For those river sections laid on landscapes where drainage is more challenging to be extracted, such as those with larger positional error, the corresponding accuracy improvements are proportionally larger.

Sensitivity analysis
The most sensitive parameter in the RSSB method is the burning intensity, φ, that controls the extent of lowering the water body and raising the land.Fig. 8a shows the probability distributions of accuracy improvement under four φ settings, ranging from 1 m to 20 m.With increasing φ, the distributions are more right-skewed, implying greater improvement in accuracy.However, when comparing average improvement (margin), increasing φ from 10 m to 20 m does not  Z.Wang et al. continue the improvement, but decrease the margin from 112.69 m to 111.18 m.Stream burning can rectify drainage location but may compromise subsequent terrain analysis (Callow et al. 2007).To preserve the observational water presence with less modification on DEM, we conducted the Pareto optimum of φ over the entire Lancang-Mekong River basin.The Pareto optimum (Fig. 8b) shows the accuracy improvement added to each previous φ setting (e.g., category φ = 5 is compared to category φ = 1, while category φ = 1 is compared to φ = 0), in relation to different DEM modification percentages.Category φ = 1 that is compared to the non-burning setting φ = 0 contributes the highest added accuracy improvement, suggesting that stream burning is vital.The optimum of φ is achieved near 10 m, considering both less DEM modification and added accuracy improvement.This value is similar to that in a previous study (Callow et al. 2007) which reported that stream burning with 10 m trenches has the least impact on subsequent terrain analysis.
To further examine the heterogeneous burning intensity, we conducted sensitivity analysis of burning intensity per stream width group.As presented in Table 2, φ = 10 and φ = 20 presents the highest accuracy improvements for all width groups.According to the sensitivity analysis (Table 2), we used the eq.6 to study the effects of heterogeneous burning intensity.Generally, it is found that the effects of heterogeneous burning intensity are not significant (Fig. 9).Both extractions with heterogeneous burning intensity (RSSB-D8 and RSSB-D∞) present marginally higher positional difference than the extraction with the fixed burning intensity (RSSB-fixed) for all river types.The resultant differences between RSSB-D8 and RSSB-D∞ extractions are also minimal.

Importance of remote sensing inputs
The methodological advance of this study mainly comes from the integration of water presence information derived from optical remote sensing.This integration enables high-resolution and multitemporal    Getirana et al. 2009;Kenny and Matthews 2005;Lehner and Grill, 2013;Lindsay 2016;Saunders and Maidment 1996;Woodrow et al., 2016;Wu et al. 2019;Yadav and Hatfield 2018;Yamazaki et al., 2019) that water map inputs are critical in order to precisely construct drainage networks in large river basins.Previous studies mostly used produced maps (e.g., blue lines on maps) to improve channel locations.Among them, MERIT-Hydro (Yamazaki et al. 2019) used water maps derived from Landsat imagery by using the conventional stream burning method, which is the only study thus far to explicitly use satellite imagery at the basin scale or larger.However, it used the mapping products rather than directly integrating satellite imagery.As a result, previous studies did not derive as much benefit from the spatiotemporal advantages (e.g., higherresolution, multitemporality) of satellite imagery as is possible for drainage networks extraction.This study, to the best of our knowledge, is the first to explore the spatiotemporal advantages of satellite imagery with a focus on the high-resolution aspect of extracting drainage networks at large scales.
In comparison with previous water maps inputs, direct integration of optical remote sensing imagery can offer many advantages.First, spatially contiguous burning delivered by direct integration facilitates the accurate extraction of drainage networks at large scales.As previous stream burning approaches often adopted blue lines from maps, which are in binary form (usually 1 for water, 0 for non-water) (Lindsay 2016;Wu et al. 2019), the thin and isolated burning can result in inaccurate connections of reaches or even disconnections (Habib et al. 2018).Our parallel experiment confirms that when a GRWL layer (river centerlines) in binary form is used as the input with equivalent extraction settings, many disconnections are produced, notably, on the main stem of Lancang-Mekong River near the dam or over low-relief areas.For drainage networks extraction over a large region, disconnections are usually inevitable due to the uncertainties in topography (Schwanghart et al. 2013;Wu et al. 2017), which could introduce more disconnections as resolution increases.However, direct integration by the remote sensing water index can generate a connected, precise drainage networks without necessitating further processing (e.g., Poggio and Soille 2012;Yamazaki et al., 2019), as demonstrated for the Lancang-Mekong River basin in this study.Its success may be attributed to spatial burning where the quantification of water presence shows a gradual decrease from river center (pure water pixel, mostly) to river edge (mixed pixel) and to bank (pure land pixel), which ensures smooth lowering (or raising for land) on the topography.Although the disconnectivity issue is not an inherent limitation of conventional methods and RSSB is not theoretically immune to disconnection (e.g., in the case of inadequate parameter setting), the spatially contiguous burning has been useful in minimizing this issue.Furthermore, the direct integration by a water index, as in RSSB, is a more objective and fully automatic method to produce drainage networks.
An additional comparison with hydrography in the selected subbasin of the Mekong River (Fig. 10) not only emphasizes the importance of remote sensing inputs but also underscores the improvements derived from spatially continuous burning and the higher resolution proposed in this study.Meanwhile, it reaffirms the importance of vertical accuracy of DEM data.MERIT-Hydro is composed of hydrography data rather than drainage networks data, which does not include some key networks information such as stream order.Hence, a more systematic comparison (e.g., river length distribution statistics as in Table 1) is not possible.Moreover, water index MuWI as the first derivative of reflectance is used in this study to identify likely waterbodies, while the second derivative of reflectance embedding river morphology can be a potential replacement, such as the singularity index which is useful for estimating river locations (Isikdogan et al. 2018(Isikdogan et al. , 2017a(Isikdogan et al. , 2017b)).
Besides, the ability to derive dynamic drainage network is determined by the multi-temporal remote sensing inputs.Multi-temporal remote sensing inputs temporalize the hydrologically-conditioned DEM within the framework of RSSB, which could transform crosssectional analysis to panel analysis of the river networks.For example, the analysis of patterns of lakes connected to river networks (Gardner et al. 2019) could possibly scale up in time dimension and to larger scales.Similarly, multi-temporal RSSB could improve hydrological modeling.For example, the concept of distributed time variant gain modeling (Xia et al. 2005) could be complemented by RSSB for its flow routing process.However, implementing time-variant flow routing may incur complications, such as the need for revising model structure and considerable additional computations.Given the existence of multitemporal DEM sources (e.g.TanDEM-X data, multi-temporal LiDAR collection), multi-temporal DEMs are not as useful as the RSSB over large scales considering their low accessibility and limitation on surface water elevation.
Moreover, the direct integration of remote sensing imagery, or the RSSB method, delivers potential for higher resolution in drainage networks extractions.As demonstrated in Section 3.1, such higher resolution is both feasible and effective for representing small rivers.Generally, higher-resolution input is not limited to Sentinel-2, or even optical remote sensing, as used in this study; other types of imagery can also be the potential sources, such as 10-m Sentinel-1 SAR (Synthetic Aperture Radar) and 3.7-m WorldView-3 SWIR (Short-wave Infrared) imagery.However, various potential challenges exist when adopting other types of imagery, such as the development of an adequate quantification on other imagery as MuWI on Sentinel-2 for optimal integration.

Linking hydromorphology with hydrology
Delineation of drainage networks can be sourced from hydromorphology, or fluvial geomorphology.Hydromorphology approaches the structure and evolution of Earth's water resources by focusing on morphological characteristics, such as the width, depth, and longitudinal profile of a river, that reflect the dynamic, historical evolution of drainage basins induced by both natural and human influences (Chen et al. 2019;Vogel 2011).The combination of hydromorphological and hydraulic characteristics, such as elevation, slope, flow direction and catchment area, has filled many knowledge gaps regarding drainage networks Z. Wang et al. properties (Allen et al. 2018;Allen and Pavelsky, 2018;Altenau et al. 2019;Barefoot et al. 2019;Frasson et al., 2019;Gardner et al. 2019;Griscom et al. 2017), particularly by efforts of the Surface Water and Ocean Topography (SWOT) mission.Further advances may be hampered by the separation of hydromorphologic measurements from drainage networks, particularly when studying them from systematic perspectives.From the hydromorphological aspect, recent progress is largely motivated by satellite remote sensing, which can provide direct multitemporal measurements at large scales.From the hydrological standpoint, flow networks topology and its quantities are still limited in terms of spatial and temporal representations.
In these regards, direct integration in this study can serve as the foundation to transform the understanding of the relationship between river networks and hydraulic quantities from that of related analysis to one of integrated fusion.Related analysis is an intuitive method that takes full advantage of current high quality drainage networks data (mostly HydroSHEDS) and relates it to satellite-derived hydromorphology, usually involving Spatial Join geoprocessing (Lin et al., 2019).The spatial discrepancy between the remotely sensed river and DEM-derived flow path is a fundamental problem in related analysis.From our calculations, the mean spatial discrepancy between HydroSHEDS and GRWL centerlines in the Lancang-Mekong River basin is greater than 700 m.This issue becomes more significant when the differences in spatial resolutions between the two datasets are larger, which poses higher uncertainties.Moreover, temporal discrepancy usually exists in related analysis because remote sensing imagery is often multitemporal, while topography is static.It is certainly appreciated that related analysis has proven feasible for many analyses (Frasson et al., 2019;Gardner et al. 2019;Lin et al., 2019), but those analyses are less sensitive to the spatial discrepancies and detailed flow information (e.g., the abundance of lakes connected to river networks (Gardner et al. 2019)), while some others are more sensitive.For example, Manning's equation-one of the key theoretical equations for deriving discharge from satellites-would be augmented at the basin scale with dense hydromorphologic measurements by combining flow propagation over a more authentic river network (Lin et al., 2019); this is important to further the understanding of intermittent and ephemeral streams in drainage networks that are currently poorly understood.
Ideally, bathymetry data is constructed over the entire basin by replacing the modeled river centerline with the real thalweg and extracting the drainage networks from the channel bed.However, this approach is still far less feasible at present due to the difficulty in penetrating water and returning signals at large scales.Nevertheless, extracting drainage networks does not necessarily require absolute elevation; it only requires relative relief between adjacent pixels.The RSSB method provides a basis to construct the quasi-bathymetry where the absolute elevation of each pixel is not guaranteed but relative relief between adjacent pixels could be improved.Given that not all wavelengths of optical remote sensing are entirely attenuated in water, reflectance could partially carry depth information under certain conditions (Alsdorf et al. 2007;Choubey 1998;Harrington et al. 1992;Legleiter et al. 2009;Legleiter and Harrison 2018).However, it is rather difficult to confirm any physical meaning of MuWI, as well as the integration in this study, which requires more comprehensive spectral analysis combined with field work.Besides, current burning is locally heterogenous.Subsequent study with the aim to reconstruct the natural basin landforms is probably achieved by a more sophisticated heterogenous burning throughout the basin.Nonetheless, our proposed method effectively extracts high-resolution drainage networks over sufficiently large areas, which may support the quasi-bathymetry assumption of the integration.In particular, its effectiveness is also revealed in its fully automated extraction without the need for postprocessing.

Limitations and future work
Despite the promising performance of our proposed integration method for the entire Lancang-Mekong River basin, there are several limitations and potential areas for future advancements in terms of both integration and flow representation.First, sophisticated remote sensing tests are needed regarding the reflectance impact, selection of burning equation, and water depth retrieval.We used Top-Of-Atmosphere (TOA) reflectance data as remote sensing inputs, as they were the highest-level reflectance data provided by the European Space Agency (ESA) over the Lancang-Mekong River basin at the time of this study (Gascon et al. 2017).However, radiometric corrections, including atmospheric, terrain and cirrus corrections which convert TOA reflectance to Bottom-Of-Atmosphere (BOA) or surface reflectance, are likely to improve inland water detection, as are water presence quantification and drainage networks extraction.A linear form of the burning equation (Eq.( 4)) was applied in this study.It assumes a linear relationship between MuWI and relative water depth, which is empirical and beneficial for the sake of brevity.A systematic examination of the extents and conditions at which water depth information can be retrieved from the remote sensing water index is still necessary in order to support the integration method.This examination will certainly contribute not only to drainage extraction but also the vertical dimension of hydromorphology, which is also essential in biogeochemical processes (Elosegi et al. 2010;Ledesma et al. 2018).
Second, potential increases in extraction resolution are possible but constrained by computation over large scales.On the one hand, the resolution of a final extraction can be further refined using the RSSB method by substituting higher resolution inputs of remote sensing (e.g., 1-m imagery over the US from the National Agriculture Imagery Program) and/or DEM (e.g., 30-m DEMs, such as SRTM30, AW3D30, NASADEM).On the other hand, it is necessary to investigate the boundary of resolution difference between the two inputs (i.e., imagery and DEM) to maintain the hydrological conditions.Furthermore, many studies on drainage extraction (Lehner et al. 2008;Yamazaki et al. 2019), including ours, identified the computational restrictions of hardware or algorithms as one of the key barriers to enhancing drainage resolution at a large scale.In the case of our 10-m resolution with 48,244 × 39,384 pixels, the computation time shows a nearly linear speedup when using parallel computing (5 nodes over a single node has a 4.7× speedup), suggesting scalability is feasible.However, we identified memory size as the major bottleneck of computation.For example, ~100 GB of memory was used to conduct the 10-m extraction of the subbasin in this study, and a theoretical >1 TB of memory would be needed for the entire Lancang-Mekong River basin.We recommend a greater number of collaborations with the community of distributed computing for the efficient parallel algorithms (e.g., Richardson et al. 2014) on more advanced high-performance computing systems.
Although our method and results have adequately integrated observational river locations, it is still insufficient in some ways to represent the physical realism of flow pathways.A braided river, in-channel flow routes, or bifurcation of rejoined or non-rejoined networks cannot be represented due to the D8 algorithm used, which is a single flow direction method.In addition, multiple flow direction methods such as D-infinity (Tarboton 1997), and recent work (Wang et al. 2020) originating from the field of graph theory may also facilitate a more realistic representation.Those multiple flow direction methods may improve positional accuracy of the extracted river locations.However, extracting river networks with explicit stream order and hydraulic connection based on multiple flow direction methods remains a challenge.For instance, to evaluate the performance of D-infinity method, we calculated flow direction and flow accumulation using D-infinity, and assessed positional accuracy for D-infinity-derived river locations (Table 3).Although the overall positional difference of Dinfinity-derived river locations is lower than that of D8-derived river networks (89.7 m vs. 123.7 m), it is recognized that upstream-downstream linkage and stream order are hard to build when using D-infinity, which is attributed to the difference between "river locations" and "river networks" in Table 3.
A 90-m DEM is used because its high vertical accuracy outweighs the effects of high horizontal resolution when compared to 30-m SRTM and AW3D DEMs.Future studies could consider DEMs with similarly high vertical accuracy of higher horizontal resolutions.As stereophotogrammetry is also useful for DEM generation (Grosse et al. 2012;Pulighe and Fava 2013), it could be feasible to explore drainage extraction solely from optical imagery in future studies.

Conclusions
In this study, we proposed an integrative method for extracting drainage networks that employs high-resolution, multitemporal optical satellite remote sensing imagery.By implementing the RSSB method using Sentinel-2 imagery over the Lancang-Mekong River basin, we demonstrated that this method is able to extract high-resolution drainage networks in contrast to the conventional model-and mapbased methods.The high-resolution extraction also exhibits improved spatial accuracy (~50% error reduction) and an enhanced ability to delineate dynamics in drainage networks.Despite some limitations, our attempt to use remote sensing to resolve spatial constraints not only contributes to drainage networks extraction methodology but also provides an improved basis for understanding the role of inland waters in global elemental cycles by refining the size and distribution of drainage network, particularly for small rivers.In addition to increased resolution, RSSB enforces the large-scale connectivity of drainage networks that is not sufficiently considered by extractions on high-resolution imagery or high-resolution DEM.The large-scale high-resolution drainage networks constructed with high accuracy could advance estimates of combined flow information (e.g., flow direction, flow accumulation, slope), networks topology (e.g., stream order, fractal, connectivity), and hydromorphological characteristics (e.g., river width, meander wavelength, sinuosity), allowing greater understanding of complex interactions between water flow and elemental exchange processes.By introducing direct integration of emerging optical remote sensing, this study presents a new method to extract high-resolution active drainage networks; this information is useful to determine flow routing in changing river systems, which is especially important for constraining estimates of smaller streams and rivers in the river network.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table 3
Comparisons of positional difference between D-infinity-derived river locations and D8-derived river networks.

Fig. 1 .
Fig. 1.(a) Lancang-Mekong River networks extracted in this study.(b) A state-of-the-art Landsat-based river centerlines extraction (GRWL) in the study area with break points (river centerline disconnections).

Fig. 2 .
Fig. 2. Schematic flowchart of the proposed RSSB method for high-resolution extraction of a drainage networks with three major steps.

Fig. 3 .
Fig. 3. Comparison of extractions in high-resolution and coarse-resolution at two zoom levels: (a) river networks level, and (b) meander level.Shown to the left of each panel are high-resolution images from Bing Aerial.Coordinates alongside images are in the unit of meter under WGS 84 / Pseudo-Mercator projection (EPSG:3857).On the right of each panel are 10-m and 90-m extraction results of flow direction, flow accumulation, and vector river line overlaid on shaded relief.

Fig. 4 .
Fig. 4. Cumulative distribution functions (CDFs) of positional difference with the observational centerlines for networks extracted by RSSB at both resolutions using various DEMs (RSSB extractions), and for the networks extracted by a conventional method (Non-RSSB extraction).

Fig. 5 .
Fig. 5. Multitemporal drainage networks (a: left panel) and water presence maps (a: right panel) at an active reach in the Lancang-Mekong River network.The temporal changes are reflected on (b) the initiation of a new river channel, and (c) seasonality between dry and wet seasons.

Fig. 6 .
Fig. 6.Improvement of spatial accuracy of extracted networks varying in the width of the river channel W. Symbol *** indicates p-value<0.01 in Kolmogorov-Smirnov test.Origin represents extractions from DEM only.RSSB represents extractions from the integration of Sentinel-2 imagery and DEM.

Fig. 7 .
Fig. 7. Spatial distribution of improvement on positional accuracy in Lancang-Mekong River basin.The improvement (D improvement , c), defined as the difference between positional errors of non-integration (D origin , a) and integration (D rssb , b), is categorized into 3 groups: Increased group for D improvement < − 10 m, Maintained group for − 10 m ≤ D improvement ≤ 10 m, and Reduced group for D improvement > 10 m.

Fig. 8 .
Fig. 8. Sensitivity analysis of burning intensity with (a) statistical distributions of accuracy improvements under four φ settings, and (b) Pareto optimum versus least DEM modification.

Fig. 9 .
Fig. 9. Effects of heterogeneous burning intensity.RSSB-D8 represents RSSB extraction with heterogeneous burning intensity calculated from D8 flow accumulation; RSSB-D∞ represents RSSB extraction with heterogeneous burning intensity calculated from D-infinity flow accumulation; RSSB-fixed represents RSSB extraction with fixed burning intensity; Origin represents the non-RSSB extraction.

Table 2
Sensitivity analysis of burning intensity per stream width group.