Coupling threshold theory and satellite image derived channel width to estimate the formative discharge of Himalayan Foreland rivers

We propose an innovative methodology to estimate the formative discharge of alluvial rivers from remote sensing images. This procedure involves automatic extraction of the width of a channel from Landsat Thematic Mapper, Landsat 8, and Sentinel-1 satellite images. We translate the channel width extracted from satellite images to discharge by using a widthdischarge regime curve established previously by us for the Himalayan Rivers. This regime curve is based on the threshold theory, a simple physical force balance that explains the first-order geometry of alluvial channels. Using this procedure, we 5 estimate the discharge of six major rivers of the Himalayan Foreland: the Brahmaputra, Chenab, Ganga, Indus, Kosi, and Teesta rivers. Except highly regulated rivers (Indus and Chenab), our estimates of the discharge from satellite images can be compared with the mean annual discharge obtained from historical records of gauging stations. We have shown that this procedure applies both to braided and single-thread rivers over a large territory. Further our methodology to estimate discharge from remote sensing images does not rely on continuous ground calibration. 10

is distributed through multiple and mobile threads (Smith et al., 1996;Ashmore and Sauks, 2006). Braided rivers are therefore often not gauged; and where these exist, the gauging stations are located at places like dams with artificially regulated flow. This hinders our ability to assess discharge in the individual threads of a braided river.
These studies establish rating relationships between some image-derived parameters (width, water level or stage, slope), to the instantaneous discharge measured in the field (Leopold and Maddock, 1953). Equations that define the hydraulic geometry of a channel relate width (W), average depth (H), and slope (S) of a channel to its discharge (Q) according to: where a, b, c, e, f, g are site specific constants and exponents. The available methods, based on remote sensing data, to estimate the discharge of a river therefore cannot be extrapolated to other rivers, or even to other locations on the same river. 35 Moreover, as these rating curves vary significantly between locations, they must be established for each location independently.
For example, Smith et al. (1995); Smith (1997); Smith and Pavelsky (2008) and Ashmore and Sauks (2006), used synthetic aperture radar and ortho-rectified aerial images to estimate discharge in braided rivers. They related the image derived effective width of a braided river to the discharge at a nearby gauge station to establish a relationship of the form of equations 1, 2, & 3.
Their approach provides an estimate of the total discharge in a braided river, at a given section. However, this technique is site 40 specific and assumes that the river bed does not change over time.
Few attempts have been made to overcome these limitations; for example Bjerklie et al. (2005) used aerial orthophotographs and SAR images to estimate discharge in various single-thread and braided rivers. To estimate the discharge they extracted the maximum water width at a given river reach. They then combined the image-derived channel widths with channel slopes obtained from topographic maps, and a statistical hydrologic model. They reported standard errors of 50 − 100 %. However, 45 after using a calibration function based on field observation, the error reduced to values as low as 10 %. Later, Sun et al. (2010) used Japan Earth Resource Satellite-1 (JERS-1) SAR images to measure the effective width of the Mekong River at the Pakse gauging station in Laos. They used rainfall-runoff model to estimate the discharge from the image-derived width and suggested that using this procedure, the discharge could be estimated in any ungauged river basin within an acceptable level of accuracy.
They established a close agreement between the measured discharge of the Mekong River at Paske station and the model 50 estimate to the 90 % uncertainty level. As discussed earlier by Bjerklie et al. (2005), later Sun et al. (2010) indicated that the precision can be improved by calibrating the rainfall-runoff model with a hydraulic geometry relation, and that a calibrated rainfall-runoff model can be used to estimate the discharge in any ungauged river using the measured width only. Gleason and Smith (2014) have suggested that the discharge of a single-thread river can be estimated from satellite images only, without any ground measurement. They plotted the exponents and coefficients of hydraulic regime equations established at 88 different 55 gauging stations along six rivers in the United States, and found that the exponents and coefficients are correlated. Recently Kebede et al. (2020) have used Landsat images to estimate daily discharge of the Lhasa River in the Tibetan Plateau. They have used image derived hydraulic variables to compute the discharge by using modified Manning equation and rating curves established from the in-situ measurement of width and discharge.
The studies discussed above attempt to address the issue of site-specificity, and propose methods to estimate discharge 60 without empirical calibration. However Bjerklie et al. (2005), and Sun et al. (2010) also show that a better accuracy in discharge prediction can only be achieved with some calibration to ground measurements. Therefore, a physically robust method to resolve the site-specificity of rating curves remains to be described.
To address this issue of site-specificity, we have developed a semi-empirical width-discharge regime relation based on the threshold theory and field measurement of various braided and meandering rivers on the Ganga and Brahmaputra plain 3 Morphology of alluvial river Lacey (1930) was the first to observe a dependency of width of an alluvial river on its discharge. Based on measurements in various single-thread alluvial rivers and canals in India and Egypt, he found that the width of a river channel scales as the square root to the discharge (e ∼ 0.5 in Eq. 1).

90
To explore the physical basis of Lacey's observation, Glover and Florey (1951) and Henderson (1963) developed a theory based on the concept of threshold channel. According to this theory, with a constant water discharge, the balance between gravity and fluid friction maintains the sediment at threshold of motion, everywhere on the bed surface. This mechanism sets the cross-section shape and size of a channel. The resulting width (W ) -discharge (Q w ) relationship in dimensionless form reads (Seizilles, 2013;Gaurav et al., 2014;Métivier et al., 2016Gaurav et al., 2017): where Q * = Q w /(d 2 s √ gd s ) is the dimensionless water discharge, d s is the grain size, ρ f ≈ 1000 kg m −3 is the density of water, ρ s ≈ 2650 kg m −3 is the density of quartz, g ≈ 9.81 m s −2 is the acceleration of gravity, C f ≈ 0.1 is the Chézy friction factor, µ ≈ 0.7 is the Coulomb's coefficient of friction, K(1/2) ≈ 1.85 is the elliptic integral of the first kind, and θ t is the threshold Shield's parameter that depends on the sediments grain size. The typical grain size of the sediments of the 100 Himalayan Foreland rivers is order of d s = 100 − 300 µm. Thus the dimensionless grain size D * = (d 3 s gρ 2 s /η 2 ) 1/3 1 − 6, where η ≈ 10 −3 Pa.s is the dynamic viscosity of water. In this range of values the threshold Shield number is on order of θ t ∼ 0.1 with a maximum around 0.3 (Julien, 1995;Selim Yalin, 1992). Recently Delorme et al. (2017), obtained an experimental value of θ t ∼ 0.25 for silica sands of size 150 µm. Here we have taken the upper value of θ t = 0.3 as a conservative estimate.
Taking lower values of threshold Shield parameter, such as the classical 0.1 would lead to a slightly better match between the theoretical prediction and the data but it does not lead to a significant change in our conclusions.
Eq. 4 is the theoretical equivalent to the Lacey's law. This theory explains the mechanism how a single-thread alluvial river, at threshold of sediment transport, adjust their geometry in response to the imposed water discharge. Strictly speaking, mean equilibrium geometry of a natural alluvial channel is not set by a single discharge, rather a range of discharges is responsible for determining the channel form (Leopold and Maddock, 1953;Wolman and Miller, 1960;Blom et al., 2017;Dunne and 110 Jerolmack, 2020). However, what value corresponds to the channel forming discharge of an alluvial river remains a matter of debate? Wolman and Miller (1960); Wolman and Leopold (1957); Phillips and Jerolmack (2016) proposed that the bankfull discharge and discharge associated with a certain frequency distribution can be used to define the channel forming discharge.
Since threshold theory predicts the morphology of a single-thread channel, one may use it to estimate discharge that relates to present day geometry of alluvial channels. To test this, we use the regime curve that we established from threshold theory 115 and measurement of hydraulic geometry of various sandy alluvial rivers in the Himalayan Foreland (Gaurav et al., 2014(Gaurav et al., , 2017. Figure 2 suggests that the individual thread of the Himalayan Foreland rivers share a common width-discharge regime relation, and to the first order their morphology can be explained by threshold theory. The theoretical exponent accords with the empirical exponent of the width-discharge curve. However, the threads are wider than predicted by a factor of about 2 ( Fig. 2). We further adjust the prefactor to the data while keeping the theoretical exponent to establish a generalised semi-120 empirical "width-discharge" regime relationship for the Himalayan Foreland rivers. This regime curve is then used to estimate the discharge of various Himalayan rivers by measuring their width from satellite images.

Dataset
To measure the width of a river channel, we use images acquired from Landsat Thematic Mapper (TM), Landsat 8 and Sentinel 125 1A satellites (Appendix A1). All images of the Landsat and Sentinel satellite missions are freely available and they can be downloaded from the US Geological Survey (https://earthexplorer.usgs.gov) and Alaska Satellite Facility (https://www.asf. alaska.edu/sentinel) websites. We have downloaded all available cloud-free Landsat satellite images, at the locations that were near the in-situ measurement stations for which discharge data was available with us ( Fig. 3). Only a few cloud free Landsat images are available for the period of June to September. This is mainly because of the strong monsoon that causes intense 130 rainfall and dense cloud covers. To overcome seasonal effect and fill the data gap during the monsoon period, we use Sentinel 1A product. Sentinel-1 satellite mission is equipped with Advanced Synthetic Aperture Radar (ASAR) sensor that operates in C-band (5.4 GHz) of microwave frequency (Schlaffer et al., 2015;Martinis et al., 2018). Advanced Synthetic Aperture Radar system can operate both day and night and has the capability to penetrate clouds and heavy rainfall. This special characteristic of SAR sensors enables uninterrupted imaging of the Earth's surface during the bad weather conditions as well.

135
In-situ measurements of average monthly discharge for some time intervals of varying length between 1949-1975 are available for the Brahmaputra, Teesta, Ganga, Chenab, and Indus rivers of the Himalayan Foreland. They can be freely downloaded

Width extraction
Our main objective is to extract the width of individual river channels from satellite images. We have developed an automated program in python 3.7 that takes a gray scale image as an input to classify the image pixels into binary water and non-water 145 classes. The pixels classified as water are the foreground object and will be used to define river channels. Dry pixels serve as a background object. To extract the river channels, we use the infra-red bands of Landsat-TM and Landsat-8 images. In Landsat-TM, the infra-red (0.76 − 0.90 µm) wavelength corresponds to band 4 whereas, in Landsat-8 image, it corresponds to band 5 (0.85 − 0.88 µm). Theoretically, since water absorbs most of the infra-red radiations it appears dark, with an associated brightness value close to 0. This typical characteristic of the infra-red signal allows a clear distinction between the water 150 covered dry areas on the satellite images (Frazier et al., 2000). However, in the case of a river, the pixel intensity varies widely because of heterogeneous reflectance of river water, due to the presence of sediment and organic particles (Nykanen et al., 1998). Because the image intensity is not exactly 0 or 1, we introduce a threshold intensity to classify the pixels. Based on this criteria, we convert the gray scale image f (x, y) into a binary image g(x, y), which distinguished between the water-covered and dry areas. This approach takes an object-background image and selects a threshold value that segments image pixels into 155 either object (1) or background (0) (Ridler and Calvard, 1978;Sezgin et al., 2004). We apply the algorithm proposed by Yanni and Horne (1994) to obtain the threshold value iteratively. Once this optimal value is obtained, we apply it to classify our pixels into water and dry classes (Fig. 4). The binary classification of satellite images into water and dry pixels can produce spurious features as well (Fig. 4). These consist of wet pixels that get classified as dry 160 or of isolated water pixels that appear randomly in the binary images (Passalacqua et al., 2013). Clusters (usually 2-3 pixels in size) that appear inside the river network do not correspond to bars or islands. We found frequent areas where strong reflection from the bed sediment cause water pixels to appear more like sand. Isolated water pixels that do not belong to the river are located in water-logged areas. We identify these types of errors and reprocess the binary images to remove them automatically.
For this, we first identify the isolated water patches from the binary images. To do this, we define a search window of 7×7 pixel 165 size. We run this window on the image and look for neighboring water pixels in all surrounding directions. If a water pixel in the classified image is disconnected in all directions from the neighboring water pixels for more than seven pixels, we consider them as isolated water bodies. We therefore re-classify such pixels as dry. We re-iterate this procedure by applying a region growing algorithm (Mehnert and Jackway, 1997;Bernander et al., 2013;Fan et al., 2005). For this we initially select a water seed pixel inside the river channel. The algorithm uses the initial water pixels and starts growing. This procedure removes all 170 isolated water patches from the binary image, and retains only water pixels connected to the river network.
Once images are reclassified, we reprocess them to merge the water pixels that were initially classified as dry inside a river channel. For this we define a search window of 3 × 3 pixels. We choose this size by assuming that dry pixels should be more than 90 meter in size to be considered as bars or islands. Otherwise, such pixels are treated as water pixels. We move the search window on the binary image and look for neighboring dry pixels inside the river channel.

175
Similarly, to identify river pixels from Sentinel 1A images, we use VH (Vertical transmission and Horizontal reception) polarized band. We have Sentinel Application Platform (SNAP) v6.0 to perform the radiometric calibration, speckle noise reduction using refined Lee filter and terrain corrections and finally generate the backscatter (σ 0 ) image. In microwave region, open and calm water bodies exhibit low backscatter values due to high specular reflection from the water surface (Schlaffer et al., 2015;Twele et al., 2016;Amitrano et al., 2018). We manually set a threshold value to separate water and dry pixels from 180 Sentinel-1 images. Finally, we follow a similar procedure as we developed for Landsat images to process the binary image obtained from Sentinel-1.
Once the satellite images are classified, we use the binary images to extract the width of each channel. We do this by measuring the distance from the center of a channel to its banks orthogonally to the flow direction. A detailed procedure of width extraction of a river channel is given in Appendix B.

Accuracy assessment
To assess the precision with which we can estimate the discharge of a thread, we need to quantify the accuracy of our widthextraction procedure using Landsat and Sentinel-1 satellite images. To evaluate this, we superimpose the contours of river channels, extracted using our algorithm, to the original gray-scale images used for the extraction. We then carefully check for 190 a match between the contours boundary and water boundary in gray scale image. We observed a good agreement between automatically extracted channel boundary and the edge of the water line in gray image. However, our algorithm fails to extract the contours of the smallest channels (60 -90 meter in width). Several reasons explain this limitation. First, as these channels are both shallow and only a few pixels wide, their pixel intensity is close to the pixel intensity of dry areas. Therefore, the optimal threshold applied to categorize the image pixels does not identify these channels as water. Second, although an increase in 195 the classification threshold could force the algorithm to identify these pixels as water, it would also add significant noise by classifying many dry pixels as water pixels. Such a limitation appears to be closely related to the image resolution.
Given this qualitative agreement, we proceed to evaluate the accuracy of the width extraction procedure. To do this, we overlay the transects used by the algorithm to measure the width of a thread on the original image (Fig. 5 a). We then manually measure the width at randomly selected transects for comparison. For each river, we manually measure the width at more than  There are some outliers however. They correspond to places along the threads where our automated procedure draws erro-205 neous transects (Fig. 5 b). Most of such transects are located near highly curved reaches at the confluence or diffluence of two or more threads. In such places, the width of a thread is overestimated sometimes by more than 50 % compared to the width measured manually. At most locations though, our procedure extracts valid transects (Fig. 5 a). Further, we assess the distribution of relative discrepancies between automatically and manually measured widths (Fig. 7).
We observe that the relative error of 90 % of our measurements is centered around a mean µ ≈ −0.02 with a standard deviation 210 σ ≈ 0.07. This validates the width-extraction procedure.

Width variability along a thread
Particularly in a braided river, the width of a thread varies significantly along its course. To quantify this variability, we select a reach and plot the probability distribution of the width measured across different transects. We observed that the distribution of width histograms is skewed Figure 8. This skewness results from the natural variability of width along the course and also due 215 to the error in width extraction from images, particularly at the location where the curvature of a thread is high. The resulting skewness will be amplified in the discharge histogram because of the non-linear relationship relating the two variables. To take the skewness into account, we have used the most probable width W m as the representative value of the width (Eq. 1). This value corresponds to the geometric mean of all measured values. However, in meandering rivers where the variability in width within a reach is not much, arithmetic mean can be considered a representative width.

Discharge estimation
We now proceed to estimate discharge (Q w ) for the Himalayan Rivers based on their channel widths measured from satellite images. To have a meaningful comparison between the image derived discharge and the corresponding in-situ measurement, we select a reach about ten times longer than the width of a river on satellite images. In the case of a braided river, we consider the widest channel to define reach length. In the selected reach we assume that discharge is conserved, there is no significant 225 addition or extraction of water in the river.
To estimate discharge of the study reach, we use a regime relation established by Gaurav et al. (2017) based on threshold channel theory (Eq. 4) and field measurements of channel's width and discharge on the Ganga-Brahmaputra Plain. The resulting regime relation is governed by:

230
where α is the best-fit coefficient, an empirical value obtained from fitting the prefactor of the regime curve (Eq.4) and W m is the most probable width. We use Eq. 6 to calculate the discharge for threads of known width. Because the river width scales non-linearly with discharge, regime relations obtained refer to the total width in the case of a braided river; and will not be the same as those obtained for individual threads. Since most of the studied rivers are braided, we first calculate the discharge for individual threads across 235 a given section. We then sum the discharge of the individual threads across a transect to compute the total discharge at a section.

Estimated Vs. measured discharge
Once monthly discharges for all the rivers are estimated from satellite images, we compare them with the average monthly discharge measured at the corresponding gauge stations. To do so, we plot the hydrographs of the estimated and measured discharges together (Fig. 9). We observe that the estimated discharge from satellite images remains constant throughout the 240 year, except during the monsoon period (June-September) when all the rivers show a significant rise. During the non-monsoon period (October-May), estimated discharges for most of our rivers are overestimated. To the first order, our approach is able to capture the rising trend of discharges during the monsoon period, however the estimated discharges are lesser than the measured discharges. Table 2 compares the estimated and measured discharges during the monsoon period. For most of our rivers, the difference between measured and estimated discharges are less than 50 %; though this difference is comparatively 245 high for the Indus (72 − 78 %) and Chenab (36 − 67 %) rivers (Table 2).

Discussion
It is important to note that the discharge estimated from satellite images does not correspond to an instantaneous discharge.
To understand the emergence of constant hydrograph from the estimated discharge derived from satellite images we explore the concept of channel forming (formative) discharge i.e. a discharge that sets the geometry of alluvial river channels. Several 250 workers Inglis and Lacey (1947); Leopold and Maddock (1953); Blench (1957), have shown that the geometry of an alluvial channel corresponds to a formative discharge (see table A3 in appendix for the definition of different discharges). They have discussed how a limited range of flows are responsible for shaping its channel. At low-flow discharge, the water simply flows through the threads without affecting their geometry. Schumm and Lichty (1965) used the concept of time span (geologic, modern and present) in defining the interrelationship between dependent and independent variables of a river system. According 255 to them, morphology of a river channel is set in the modern time span (last 1000 years) by the average discharge of water and sediment. In the present time span (1 year or less), channel morphology can be considered as independent variable against instantaneous discharge of water and sediment.
Similarly, it has been argued by Inglis and Lacey (1947); Leopold and Maddock (1953); Blench (1957) that it is not the highest flows that contribute the most in shaping a river channel. Such high discharges are capable of transforming the channel, 260 but they occur so infrequently that, on average, their morphological impact is small. Wolman and Miller (1960) highlighted that the bankfull discharge that occurs once each year or every two years sets the pattern and channel width of the alluvial rivers.
Formative discharge for the Himalayan rivers is expected to occur in the monsoon period, thus one may expect that during low flow such rivers maintain their flows without modifying the existing channel geometry (Roy and Sinha, 2014). This clearly reflects in the discharge hydrographs estimated from the measurement of channel's width from the satellite images ( Fig. 9.

265
Furthermore,  have recently shown that non cohesive streams laden with sediments cannot have a width much larger than the width of a threshold stream before they start to braid. They also showed that, for experimental braided rivers, threads are always formed at the bankfull flow, and at the limit of stability. Our hypothesis is thus that the formative discharge of threads in the Ganga plain is the bankfull discharge." This is probably why our estimated discharge from satellite images remain constant throughout and is mostly overestimated than the measured discharge at gauge stations.

270
According to Inglis and Lacey (1947), rivers approach their equilibrium geometry for a formative discharge that approximately corresponds to the bankfull discharge. They suggested this discharge lies between 1/2 and 2/3 of the maximum discharge. It has also been suggested that the formative discharge corresponds to the median discharge (Blench, 1957). In their study, Leopold and Maddock (1953) used the discharge that corresponds to a given frequency of occurrence and compared it to the hydraulic geometry of the river. Based on their observations in the United States, they recommended the use of the annual 275 average discharge as a proxy for the formative discharge. Hereafter, we use the definition of Leopold and Maddock (1953).
Based on our understanding of the geometry of alluvial river channels, we argue that the width of the thread that we extract from satellite images corresponds to a formative discharge. Therefore for a given river, discharge estimated from these widths should compare with the formative discharge. We now evaluate how the discharge estimated from satellite data varies with time.
We plot the monthly discharge estimated for all of our rivers to their corresponding average monthly discharge measured at the 280 gauge stations (Fig. 9). The monthly average discharge of the Himalayan Foreland rivers appears to be a representative of the actual hydrograph (Fig. C1). As suggested earlier by Inglis and Lacey (1947); Leopold and Maddock (1953) and Blench (1957) we observe that the estimated discharges from images are nearly constant throughout the year, with only small fluctuations around their mean. This supports the hypothesis that the width of the thread extracted from satellite images corresponds to a formative discharge.

285
Now we compare the discharge estimated from satellite images to the discharge measured at a nearby gauging station. To do this, we first compare the annual average discharge estimated from Landsat and Sentinel-1A images for different months to the annual average discharge measured at corresponding ground stations. We plot these discharges on a log-log scale (Fig. 10).
The discharge estimated from satellite images agrees to an order of magnitude with the measured discharge. The difference between measured and estimated annual average discharges for the Brahmaputra,Ganga,Kosi,and Teesta 290 rivers is less than 20 %. However, this difference is comparatively high for the Indus (78 %) and Chenab (49 %) rivers. Interestingly, the estimated discharge for the Teesta (at Anderson station), Ganga (at Farakka & Paksay) and Brahmaputra (at Bahadurabad) rivers converge to their measured discharge with a small difference of 5 %, 8 %, 4 % and < 1 %, respectively (table 2); whereas the estimated discharge of the Teesta (at Kaunia station) and Kosi (at Bhimnagar) show a relatively higher difference of 19 %, and 16 %. This difference could be possibly related to the anthropogenic impact on the natural flow condi-295 tion. For example, the selected study reaches for the Teesta (at Kaunia station) and Kosi (at Bhimnagar) rivers is located near the barrage where flow is highly regulated. However, this relationship is not entirely clear at this stage.
Similarly, the large difference between the estimated and measured discharge of the Indus and Chenab rivers could possibly be related to a series of dams and barrages (Kotri barrage, 1955, Mangla dam, 1967, Tarbela dam, 1976) that have been constructed. Such interventions have significantly altered the water and sediment discharge of the Indus river. For example, 300 downstream of the Kotri barrage, the average annual water discharge of the Indus river has declined at an alarming rate of about (107 to 10) ×10 9 m 3 from 1954 to 2003 (Inam et al., 2007). This continuous decline in the average annual discharge might have significantly modified the geometry of the Indus river. Further to understand our estimates of discharge for the Chenab, Indus and Teesta rivers, we plot their monthly discharge time series recorded at the corresponding gauge station together with the discharge estimated from satellite images ( Fig. C2 and 305 C3). Despite a large variability, the discharge time series of Indus and Chenab rivers show a strong declining trend during the monsoon period (June-September); whereas discharge during the non-monsoon period appears to remain constant around the mean value. Figure (C2) clearly show that discharge estimated from satellite images plot within the variability of the observed trend. The estimated discharge of the Teesta River also plot within the noise of the observed trend. This gives a confidence in our estimates of discharge, especially for the rivers we have the limited and old record (1973)(1974)(1975)(1976)(1977)(1978)(1979) of in-situ discharges.

310
In a recent study Allen and Pavelsky (2018) measured the width of the global rivers from Landsat images for the month when they commonly flow near mean discharge. We have used the water mask binary images from "Global River Width from Landsat (GRWL) database" and measure the threads width of the Brahmaputra, Chenab, Ganga, Indus and Teesta rivers (Appendix C). We used these widths to estimate discharges using our regime curve and compare them with the mean annual discharge recorded at the corresponding gauge station and our estimates from satellite images. We observed for most of our 315 rivers the discharge estimated from thread's width extracted from the GRWL database of Allen and Pavelsky (2018) falls within the same order of magnitude to the yearly average discharge measured at the corresponding gauge stations (Table C1).
This suggests that water mask from GRWL database can be used as a first order approximation of the mean discharge of the Himalayan Foreland rivers. Also we noticed that the discharge estimated from GRWL database appears to occur during early (June, July) or post (September, October) monsoon period. can be obtained easily from field measurements in the field. Since our regime equation is established from measurement of a wide range of channels spanning over the Ganga and Brahmaputra plains, we believe it can be used to obtain the first order estimate of formative discharge of rivers in the Himalayan Foreland by just measuring their channel width on satellite or aerial images.
Using our semi-empirical regime equation and satellite images of Landsat and sentinel-1 missions, we have estimated the 335 discharge of six major rivers in the Himalayan Foreland (Brahmaputra, Chenab, Ganga, Indus, Kosi, and Teesta). Our estimated discharges closely compare with the average annual discharge measured at the nearest gauging stations. This first-order agreement although encouraging, requires further research to improve the degree of agreement between measured and estimated discharges. One of the main source of uncertainty in discharge estimate is due to the error in the measurement of thread's width. This depends on the image resolution and accuracy of the algorithm used for extraction of river pixels from remote 340 sensing images. A better resolution remote sensing images would most likely minimise the uncertainty and improve the agreement between estimated and in-situ discharge. Further our regime equation established for Himalayan rivers is based on a simple physical mechanism that explains the geometry of alluvial channels. We therefore suspect that the procedure we have established could be extended to most alluvial rivers. Globally it has been observed that the threshold theory well predicts the exponent of the regime equation (Eq. 6), however the prefactor may vary significantly depending on the grain size distribution, 345 turbulent friction coefficient and the critical shield parameter . It is therefore suggested to modify this regime curve from the measurement of width, discharge and grain size of a individual threads of alluvial channels in the field before applying it to the rivers of different climatic regime. Further it should be noted that our regime curve relates to the measurement of hydraulic geometry of individual threads of braided and meandering rivers, therefore it is applicable only at the thread scale. Since the resulting regime curve is non linear, estimating discharge across a transect in a braided river from 350 the aggregated width will be different from the one obtained after the summation of discharges of the individual threads.
This study presents a robust methodology and is a step towards obtaining first order estimates of formative discharge in ungauged river basins solely from remote sensing images. It can be used for the sustainable river development and management to ensure regional water security, especially in the regions where river discharge is not readily available.

A1 Satellite images
Detailed specification of satellite data (Landsat and Sentinel-1) used in this study.

B0.2 Removal of artefacts
Thresholding a gray scale input satellite image into binary class (water and dry pixels) produces spurious features. These consist of wet pixels that get classified as dry or of isolated water pixels that appear randomly in the binary images. Clusters

375
In-situ data (discharge)  Table A2. Satellite images used for the extraction of channels width. In-situ discharge data is freely available and were downloaded from (http://www.rivdis.sr.unh.edu/maps/asi/).
(usually 2-3 pixels in size) that appear inside the river network do not corresponds to bars or islands. They appear to be more frequent in the areas where strong reflection from the bed sediment cause water pixels to appear more like a sand. Isolated water pixels that do not belong to the river are disconnected and located in water-logged areas. We have identified both of these type errors from binary image and reprocess to remove them automatically. While doing this we first define a seed point inside the main channel and run the flood filling algorithm (Mehnert and Jackway, 1997;Bernander et al., 2013;Fan et al., 380 2005). This identify water pixels in a river channel those are connected and remove the isolated water pixels those have poor connectivity (Fig. B1).

B0.3 Extraction of channel's skeleton and contour
Our channel width extraction algorithm requires to river's centerline and boundary. A river centerline often called skeleton in computer vision corresponds to its median axis. To identify the river skeleton, we have used a thinning algorithm to extract 385 river's centerline. The algorithm iteratively reduces the boundary pixels in a way that preserves its topology (for example eroding pixels must not alter the geometric properties of the object studied) and connectivity (Fig. B2 a). The final skeleton is centered within the object and reflects its geometrical properties (Zhang and Suen, 1984;Baruch, 1988;Lam et al., 1992;Chatbri et al., 2015). The thinning algorithms produces several small centre line segments, often less than 300 meter in length, that are disconnected from the channel network at one end. These segments of the skeleton are too small to be considered as 390 part of the river network. For our purpose we consider such segments as noise and filter them out. We do this iteratively, by looking for skeleton segments those are disconnected from the skeleton network at one end. To extract the channel banks, we have applied a contour extraction algorithm that detects the outer boundary of a channel (Fig. B2). The algorithm relies on a pixel-neighbourhood analysis, where a pixel in a binary image is considered a contour pixel, if it has at least one background neighbour (Chatbri et al., 2015).

B0.4 Channel's width calculation
Once the satellite images are processed to extract skeleton and channel's banks, we then proceed to extract the width of each channels. We do this by measuring the distance from the centre of a channel to its banks orthogonally to the flow. From the skeleton of the image we draw a perpendicular line to the river bank and measure the the euclidean distance (Fig. B3). In case of a braided river, especially near the junctions where more than two river join or bifurcate form a complex network. At such 400 locations our algorithm fails to measure correct width. To circumvent this we identify all the junctions from river skeleton ( Fig. B2b). In the proximity of junction we consider an area of 5 pixels and define them as a zone of channel's confluence and diffluence. In this zone we avoid to calculate the width of the channels.
Finally, we draw perpendicular transects from each pixel of the skeleton to both side of the channel and calculate the distance from any point (x, y) on the skeleton to its corresponding left (x 1 , y 1 ) and right (x 2 , y 2 ) points on the channel boundary ( Fig.   405 B3). We then sum these widths to get the total width across a transect. For simplicity, at every one kilometer distance along the channel we compute the most probable width of each channels across a river section. Finally, the discharges through a section can be calculated along an entire reach (Fig. B3). Circle in blue is the discharge estimated from satellite images.

C3 Comparison of mean annual discharge with GRWL database
Allen and Pavelsky (2018) measured the width of the global rivers from Landsat images for the month when they commonly flow near mean discharge. In their database, Global River Width from Landsat (GRWL), for braided river they have reported the aggregated width of all the active threads. This width can not be used to estimate discharge from our regime curve that 415 we established for the Himalayan Rivers. Our regime curve relates to the measurement of hydraulic geometry of individual threads of braided and meandering rivers (Gaurav et al., 2014(Gaurav et al., , 2017, therefore it is applicable only at the thread scale. Since the resulting regime curve is non linear, estimating discharge across a transect in a braided river from the aggregated width will be different from the discharge obtained from the summation of discharge of the individual threads. To overcome this, we have used binary water mask images from GRWL database to extract width of the individual threads. 420 We then use these threads to estimate their discharge using our regime curve (equations. 4 and 6 in the manuscript). We observed for most of our rivers, discharge estimated from threads width extracted from the GRWL database falls within the same order of magnitude to the yearly average discharge measured at the corresponding gauge stations (Table C1). Brahmaputra Bahadurabad 21751 ± 2942 21717 ± 4740 11149 ± 5122 Table C1. Annual average discharge measured at the gauge station and estimated from satellite images. QGRWL is the discharge estimated from binary water mask from GRWL database from Allen and Pavelsky (2018).
Author contributions. KG, AVS and FM have conceptualised the study, KG collected the field data, AVS and KG developed the algorithm to process satellite images, AK has processed the Sentinel-1 satellite images. KG has written the manuscript and FM, RS and SKT have reviewed. All authors discussed the results and contributed to the final manuscript.
Frazier, P. S., Page, K. J., et al.: Water body detection and delineation with Landsat TM data, Photogrammetric Engineering and Remote