dataesr/scientific-paragraphs-categorization huggingface dataset:
We assume that on a monthly scale, the effect of water routing is negligible for small-to mid-size basins. We sought to develop a large database of global gaged basins that would represent a range of geographic locations, environments, and basin sizes to better sample the global water cycle on Earth. We selected gages for our analysis based on data quality, geographic and temporal coverage, and location. We considered gages with an upstream area ≥2,500 km 2. We obtained river dischage data from 3 sources. First, we selected 1,737 gages from the Global Runoff Data Center (GRDC) and supplemented it with information from two other sources to fill in blank spaces on the map (notably Asia and Australia). The GRDC database contains historical mean daily and monthly discharge data from 159 countries  (WMO, 1989; BfG, 2020).
The GRDC database contained 10,361 stations when we acquired data, however the majority of these did not fit our criteria for spatial and temporal coverage. Second, we obtained data for 272 gages from the Global Streamflow Indices and Metadata (GSIM) archive  (Do et al. , 2018; Gudmundsson et al. , 2018). Finally, we obtained runoff data for 47 gages in Australia from their Bureau of Meteorology (BOM)'s Hydrologic Reference Stations (Australia BOM, 2020). We calculated monthly average runoff for months with at least 25 days of data. Volumetric flow rates in m 3 /s were converted to area-normalized fluxes in mm/mo by dividing by the land surface area in km 2 and multiplying by an appropriate conversion factor. The spatial coverage of our final 2,056 river gages (and their basins) is uneven across the globe (see Figure  1 ).
North America is over-represented with 1,111 gages (more than half the total), as is Europe with 393 gages, while we have only 70 gages in Africa, 178 in Asia, and 195 in South America. Runoff observations are limited, as they are only available at gaged locations. As an alternative, indirect estimates of runoff are available from several sources. For our experiments in closing the HC, we used estimated runoff from GRUN Ensemble  (Ghiggi et al. , 2021). The authors created a global gridded dataset of runoff with a random forest model using P and near-surface temperature as predictor variables. For the 2021 GRUN Ensemble project, the authors used input data from 21 different sources, "including a set of atmospheric reanalysis, post-processed reanalysis and interpolated-stations data. " In order to check the quality of the GRUN dataset, we performed an independent evaluation against our 2,056 gages and found that GRUN is a relatively good fit to observed discharge.
Environmental Indices We also collected observations of ancillary environmental data as inputs to our NN model. Our hypothesis is that errors in EO data (that the NN will attempt to correct) are the consequence of certain environmental conditions. For example, precipitation estimates are often biased in mountainous regions, or in relation to snow cover. The environmental data are listed with their source/citation below. In all cases, we rescaled and reprojected datasets as necessary and calculated spatial means for river basins as described above. The 12 environmental indices are: 1. Aridity index  (Trabucco & Zomer, 2019)  2. Mean elevation  (Amatulli et al. , 2018)  3. Median slope  (Amatulli et al. , 2018)  manuscript submitted to Water Resources Research 4. Basin centroid latitude (calculated) 5.
Enhanced vegetation index  (Didan, 2015)  6. Vegetation growth/senescence (calculated) 7. Irrigated area  (Siebert et al. , 2015)  8. Fire: burned area  (Giglio et al. , 2020)  9. Snow cover  (Hall & Riggs, 2021)  10. Solar radiation  (Hogan, 2015)  11. Land surface temperature  (Wan et al. , 2021)  2. 6 Preliminary Analysis of EO Datasets  and the 90%-ile). The top boxplot in each set of observations is for all pixels over land, while the lower box shows the distribution across our 2,056 gaged basins. For most variables, the distribution of fluxes is greater over the pixels compared to the basins, with higher highs and lower lows. This is particularly the case for P , but is also seen with E. manuscript submitted to Water Resources Research When we calculate the mean flux over a basin, it tends to smooth out the extremes and compress the distribution of observed fluxes.
There are also differences in the distributions within each category of fluxes. For example, GPM-Imerg contains higher observations of P , with a higher 75-and 90-percentile than the other two datasets. The monthly water storage change, ∆S, is centered at about zero for each dataset. This is expected, as the storage in pixels and basins tends to fluctuate seasonally, and any long-term trend is small compared to the annual variations in storage. Runoff has the smallest magnitude of any of the hydrologic fluxes, with a low of 0 mm/mo (no observed flow) to a 90%ile of 68 mm/mo, lower than the 90%-ile of P or E. We conclude that the lack of consensus among datasets is further evidence of the need for them to be reconciled.
Training database at basin scale over the world We analyzed the water balance at the scale of two geographic units: river basins and pixels. For the pixel-scale analysis, we used a 0. 5°equirectangular grid. We excluded Antarctica, Greenland, and the Arctic north of 77°. Near the equator, a single pixel is about 56 km on a side and has an area of about 3,100 km 2. The disadvantage to studying the water balance at the pixel scale is the lack of observations of horizontal inflow or outflow. In this regard, there are advantages to working at the scale of river basins, or watersheds. A watershed is defined as the area on the Earth's surface where water drains to a common outlet and is determined by the topography of the land surface. We obtained basin geodata in shapefile format from the GRDC, which covered many of our gaged basins  (Lehner et al.
, 2008). However, some of these basin boundaries appeared to be inaccurate. Therefore, we created a new set of boundaries for every watershed using the best-available global hydrographic data. We used a hybrid method that uses both vector-and raster-based data  (Heberger, 2022). Our method uses the vector dataset MERIT-Basins  (Lin et al. , 2019 (Lin et al. , , 2021)) , where rivers are encoded as polylines and catchment boundaries as polygons. Our resulting 2,056 basins vary in size from 2,500 km 2 to 4. 7 million km 2 for the Amazon basin. The distribution of basin sizes is highly skewed, meaning that we have manuscript submitted to Water Resources Research many small-and medium-sized basins, and fewer large basins. Of these, 119 watersheds have an area of less than 3,000 km 2. Our set of gaged basins covers 47 million km 2 , about 35% of our study domain, the land surface below 77°North and excluding Greenland and Antarctica.
For each basin, we calculated the average for each of EO variables in Tables 1 and the environmental variables described above. To calculate the spatial weighted mean, we converted each basin polygon to a grid "mask," where each pixel is a floating-point number representing the fraction of the pixel's area that is inside the basin (from 0 to 1). Because the surface area of pixels varies by latitude, we used the pixel's area in our calculation of the weighted mean. When working with a gridded runoff dataset, we are free from the constraint of using gaged basins, and we may define river basins of any size and at any location for better sampling of many environmental conditions. To take full advantage of this, we created a set of 1,698 synthetic river basins, shown in Figure  1 (b). These represent real, physical basins, but their outlets do not correspond to a gage.
We created the synthetic basins by using a gridded dataset of flow direction created by the developers of the GloFaS-LISFLOOD model  (Harrigan et al. , 2020)  and the Python library pysheds  (Bartos et al. , 2023). The basins range in size from 20,000 to 50,000 km 2 ; the relatively small size allows fits the hypothesis that we may neglect water travel time at the monthly time scale. The color coding in Figure  1  shows the experimental partition we created for the training and validation of the NN model, with 80% of basins for training (in blue), and 20% of basins for validation (in red). The result is a set of 1,698 basins shown in Figure  1 (b).
Optimal Interpolation We follow previous studies  (Aires, 2014; Pellet, Aires, Munier, et al. , 2019)  in using OI to integrate EO datasets and balance the water budget at the river basin scale. We refer to the water balance residual as the imbalance, calculated by: I = P -E -∆S -R. (3) The OI approach is based on forcing I to equal zero, or minimizing I, while distributing the errors among the inputs in inverse proportion to each variable's uncertainty. These methods, well described by  Rodgers (2000) , are referred to as "inverse methods" and are widely used in remote sensing. The goal of OI is to combine these multiple estimates to obtain the best consensus of the HC state. For a detailed explanation of the mathematics behind OI, see  Aires (2014)  and  Pellet, Aires, Munier, et al. (2019).
We begin by defining an initial best guess for each of the four HC components. This is done by calculating a weighted average of the inputs for each component. Simple weighting (SW), as described by  Aires (2014) , is an application of the method of inverse variance weighting to the problem of calculating the "best estimate" that combines multiple satellite observed fluxes. It is a form of weighted averaging where the weight on the each observation is the inverse of the variance or uncertainty of that observation. The uncertainty must be estimated a priori for each observation or dataset. After calculating the best first guess of the water budget based on the SW mean, we apply a post-filter to enforce the water balance. The post filter is a linear transformation based on the uncertainty in each component.
More generally, a NN is a flexible model that can simulate complex nonlinear relationships. Given the correct model form and proper training, it can fit any arbitrary function. Often, classical NN architecture is fully connected, meaning that every neuron has a connection with all the neurons of the previous layer. This is not the case here, where we are operating multiple independent NNs for calibration and mixture. We experimented with a number of NN architectures. While the one shown in Figure  4  is among the simpler models that we tried, it performed the best. On the left are the model inputs, the uncorrected EO datasets, and on the right are the targets, the solution from OI that results in a balanced water budget. We chose a modular architecture with separate calibration and mixture steps that allows us to investigate the outputs of individual layers as we may gain useful information from each: • First, a set of NNs serves to calibrate the individual inputs, or to transform them such that they more closely match the OI solution that satisfies the water balance constraint.
We also performed the equivalent smoothing on validation datasets in order to ensure a fair comparison. -A E 1 E 1,cal Gleam-B E 2 E 2,cal E mix ERA5 E 3 E 3,cal Total Water Storage Change CSR ΔS 1 C7 ΔS 1,cal GSFC ΔS 2 C8 ΔS 1,cal ΔS mix JPL ΔS 3 C9 ΔS 1,cal Runoff GRUN R C10 R ,.
Results Here, we evaluate the results of our optimization procedure for EO data using OI and NN modeling. The best model will be one that reconciles the inputs and results in a lower water budget residual, I. It should yield results that are plausible while changing the inputs as little as necessary. Figure  5  is an example showing the inputs and outputs of our method over one river basin. The data is for the White River at Petersburg, Indiana, United States, with a drainage area of 29,000 km 2. While no river basin is typical, this location does a good job demonstrating the output from our calculations as it has a long record of river discharge. The corrections made in this basin are relatively modest; over this region of the eastern United States, remote sensing datasets tend to be more reliable and well-calibrated due to the density and availability of in situ calibration data.
The time series plots in Figure  5  show the inputs (EO datasets, in gray), the outputs of OI (green) and the outputs of the mixture NN (purple). There is significant disagreement among the 3 P datasets as their seasonality differs. E for this location is more The NN optimization of P , E and R results in a significant improvement in I. One of the key features of our model is that it should make minimal modifications to the inputs while moving closer to a solution that balances the water budget. In particular, we note that discharge R is changed less by the NN optimization than it is by OI. This means that the NN optimization acts mostly on P and E towards a better coherency with R and ∆S. mixture model has a slight positive impact (I = -0. 03 ± 24 mm/mo). It appears therefore that most of the improvement comes from the initial calibration layer with an additional but minor improvement from the mixture layer.
Reconstruction of Total Water Storage Change An advantage to the NN architecture described in Figure  4  is that it is modular. Each step (calibration, mixture) results in an improvement to EO datasets, in terms of producing a balanced water budget, as seen in Figure  6. This is very valuable when faced with missing data: A missing HC component can be estimated by inference from the other three. Indeed, several studies have exploited this relationship (see e. g. ,  Rodell et al. , 2011; Munier et al. , 2014; Liu et al. , 2020; Pellet, Aires, Papa, et al. , 2019; Lehmann et al. , 2022). We used this approach for indirect estimation of ∆S. This allows us to estimate GRACE-like TWSC for the time period before 2002 when the satellites were launched, or to fill in missing data. If we assume a balanced water budget, rearranging Equation  1 gives ∆S = P -E -R.
Estimating missing components using indirect observations should be improved when using the optimized water components of the previous section. This is therefore an indirect evaluation of the water budget obtained by our integration framework. Overall, we obtained a significantly improved fit to GRACE observations with ∆S obtained from the three other NN-calibrated fluxes, compared to similar estimation with uncorrected EO data. At the pixel scale, our new NN-inferred ∆S compare favorably to those predicted by  Zhang et al. (2018). ment of the water components, which should be beneficial for future applications. The fact that our NN model performs well under most conditions is encouraging. Before 1980, there are fewer combinations, as some datasets have limited temporal coverage (see Table  1 ).
The green line shows TWSC from GRACE, where available (average of the three solutions in Table  1 ). The orange line is our reconstruction of ∆S. Finally, the dashed purple line is ∆S from the study by  Zhang et al. (2018)  There are also other limitations to the reconstructed datasets of TWSC. It can be shown that even a very small bias makes it impossible to calculate the trend in TWS with any degree of accuracy. We are computing TWSC from climate data only, while it has been shown that human activities like groundwater pumping and the filling and draining of reservoirs have a major impact on TWS  (Rodell et al. , 2018). We explored novel methods of analyzing and combining earth observation datasets describing major hydrologic fluxes, with the goal of reducing the overall error in estimating the water budget.
We applied a closed-form analytical solution, optimal interpolation, which forces the water budget residual to zero. This approach has several advan- The NN framework introduced by Aires (2014) and expanded upon in this paper opens new doors for the integration of satellite data to study the HC. The NN model we developed is original in the field of water budget closure studies, and has some special features that allow us to integrate satellite observations. Our model is nested, featuring independent calibration and mixture models to stay closer to the physical treatment that we intend to produce. Our approach optimizes EO datasets and closes the HC without the use of a simulation model. Rather, our data-driven approach can be set up to rely only on data from satellite returns. This makes it valuable for the calibration and validation of climate models and hydrologic models, among other applications.
consistent. The three GRACE datasets for TWSC or ∆S are highly correlated with one another, as expected since they are derived from the same satellite data. The bottom plot shows the HC residual or imbalance, I. The gray lines show each of the 27 possible combinations of the datasets (3P ×3E ×3∆S ×1R). The imbalances based on uncorrected EO data are significant: the seasonal I can reach ±50 mm/mo depending on the combination of datasets. The objective of our integration technique is to reduce this imbalance as much as possible. I based on the OI solution (in blue) is equal to zero by definition.
Over the selected basins, the reconstructed time series of TWSC do a good job recreating the seasonal patterns observed by GRACE over river basins of a range of sizes. Further, both reconstructed time series of TWSC are a significantly better fit to observations compared to estimates based on uncorrected EO data. As shown in 8, the reconstruction based on this study's NN is a slightly better fit to observations compared to the results from Zhang et al. (2018). This study's indirect estimates of TWSC are able to cover a longer time period; the modular nature of the calibration NN model allows us to use whichever dataset(s) are available in a given time period for estimation. In general, estimates are more robust when more datasets are available. As fewer datasets are available from 1980 to 2000, this is an additional source of uncertainty for hindcast estimates of TWSC.
Yokoya et al. [31]  studied flash floods and debris flows by combining (i) a numerical simulation to generate large training datasets and (ii) deep learning with Convolutional Neural Networks (CNNs) to link the synthetic binarized information to actual change images from bi-temporal remote sensing data. Most of these works, however, only dealt with single sites, impeding wider practical use. Cerbelaud et al. [27, 32]  developed the Sentinel Plot-based Change Detection (SPCD) method for exhaustive, automatic and repeatable identification of land surface damages caused by PFs in the days following an extreme weather event. Featuring limited false-positive rates (~10%), SPCD consists of a Gaussian Process (GP) object-based classifier and works very quickly at large scales (thousands of km 2 ) even on a standard local computer.
For this, several approaches were tested (unit normal global normalization, by-patch normalization). However, flattening the data was not a viable option as spectral variability was key for detecting land-cover irregularities caused by PFs. Global histogram matching was eventually the only conclusive method (not shown here) through the "match_histograms" function from the exposure module of the scikit-image Python package. It consists of searching (for each color level of a reference image channel) the level in the target image channel for which cumulative distribution functions of the two images equalize. This technique is particularly helpful to balance detector responses depending on illumination or calibration techniques.
In order to evaluate the capacity of ML/DL classifiers to be trained on Pléiades data and later function with other types of VHR sensors, orthophotos of the Vésubie and Roya valleys were upsampled to 0. 5 m spatial resolution and histogram matched with Aude 1 Pléiades' final processed image. In particular, the three red/green/blue channels of the orthoimages were converted from [0:255] to the targeted reflectance range through simple linear transformation. Emphasis was particularly placed on assessing whether the FuSVIPR method could be fed with various kinds of VHR visible data though they may originate from different sensors with contrasting wavelength ranges, optics structures, calibrations and radiometric processing techniques. Because most VHR airborne instruments do not include an NIR channel, the NIR band of Pléiades was left out in the training process to ensure transferability to other VHR products (loss in accuracy was assessed beforehand).
Multispectral Optical Change Images Furthermore, just as in Cerbelaud et al. [27] , precise change images were generated for all 4 study events from Sentinel-2's closest pre-event and post-event satellite data (Table  1 ). S-2 change information was vital to confirm actual recent modifications induced via intense surface runoff on post event VHR data. Corresponding S-2 tiles for each event were T31TDH for Aude 1 and Aude 2, T32TLP for Alpes-Mar and T36JTN for Durban. Satellite products were retrieved from the French THEIA platform for Aude 1, Aude 2 and Alpes-Mar with level-2A processing  [54]. For Durban, level-2A products were later downloaded directly from the Copernicus Open Access Hub. All of them consisted in orthorectified bottom-of-atmosphere (BOA, i. e. , corrected for atmospheric effects) flat reflectance data with sub-pixel multispectral registration.
They were acquired either from S-2A or S-2B satellites (see Table  1 ). Pre and post event images were co-registered with post event images as master images using the ENVI ® software. Each resulting pair of S-2 acquisitions presented values of average sum-of-squared intensity differences (SSD) lower than 0. 005 on all four study events. Change images simply consisted in pixel-by-pixel relative difference of vegetation-specific indices based on red (B4), green (B3) and NIR (B8) bands of S-2 (see the formula in item 5 in the left boxes of Figure  1 ). Indeed, field observations have shown that intense overland water flow is essentially characterized by the presence of eroded soils, mud or sediment deposits and deteriorations in the vegetation cover. It can thus be tracked down through vegetation-and water-based spectral variations in time and space.
In particular, NDVI  [55]  and NDWI  [56]  have proven very relevant in previous works to detect PF-induced deteriorations  [27]. The change images will be hereafter referenced to as (RD k ) S-2 , with k representing the corresponding index. With a 12 m location accuracy and a medium spatial resolution of 10 m for visible and NIR channels (Table  1 ), S-2 s revisit frequency (~5 days) was fundamental to characterize recent land cover changes that could appear ambiguous on post event images. After cubic resampling from 20 m to 10 m, the B11 SWIR band of S-2 was also considered for identification of post PF traces  [27, 57]. Indeed, it has been shown to be related to soil moisture  [58] , which is indicative of floods and a known driver of landslide occurrence  [59]. Additionally, in order to fuse VHR post event data with S-2 change images for use as predictors in pixel-based classifiers, co-registration between S-2 post event images and VHR images was also performed using the GeFolki software  [60]  prior to the generation of the S-2 change images.
All S-2 final processed change images were also histogram matched with Aude 1 S-2 data for reference using the same Python exposure module as was used for VHR images. (RD k ) S-2 values (k = NDVI, NDWI and SWIR) were finally re-mapped to the same reflectance interval as Pléiades ([0,1]*10,000 for integer storage) to avoid feeding the U-net with negative values and in order for it to account for all predictors (VHR and S-2 changes) in a rather balanced way. All products were used (and projected to, if necessary) in the native reference coordinate system of Sentinel images, World Geodetic System (WGS) 84, either Universal Transverse Mercator (UTM) zone 31N for Aude 1 and Aude 2, UTM zone 32N for Alpes-Mar and UTM zone 36S for Durban.