Xiucheng Yang

Qiming Qin

Hervé Yésou

Thomas Ledauphin

Mathieu Koehl

Pierre Grussenmeyer

Zhe Zhu

Monthly estimation of the surface water extent in France at a 10-m resolution using Sentinel-2 data

Keywords: dynamic mapping, Google Earth Engine, Sentinel-2, water bodies, France, superpixel, spectral indices, monthly

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction 1.1. Background Water bodies are the main component of the land surface. The accurate spatial detection and dynamic monitoring of inland water bodies are important tasks in many applications, such as sustainable land and water management [START_REF] Zou | Continued decrease of open surface water body area in Oklahoma during 1984-2015[END_REF], water volume and water level estimation [START_REF] Crétaux | Lake Volume Monitoring from Space[END_REF] high revisit frequency (approximately five days) using a bi-satellite system, which is important for land cover dynamic mapping and updating. Sentinel-2 can effectively extract the distribution of land surface water at a resolution of 10 m [START_REF] Du | Water bodies' mapping from Sentinel-2 imagery with Modified Normalized Difference Water Index at 10-m spatial resolution produced by sharpening the swir band[END_REF][START_REF] Gong | Stable classification with limited sample: transferring a 30-m resolution sample set collected in 2015 to mapping 10-m resolution global land cover in[END_REF][START_REF] Wang | Multi-SpectralWater Index (MuWI): A Native 10-m Multi-SpectralWater Index for accuratewater mapping on sentinel-2[END_REF][START_REF] Yang | Urban surface water body detection with suppressed built-up noise based on water indices from Sentinel-2 MSI imagery[END_REF]Yang et al., , 2017;;[START_REF] Yésou | Synergy of Sentinel-1 and Sentinel-2 imagery for wetland monitoring information extraction from continuous flow of sentinel images applied to water bodies and vegetation mapping and monitoring[END_REF], but its potential for dynamic water body monitoring and change detection deserves more attention considering the high revisit frequency of the satellite (Ogilvie et al., 2018a).

Related works

(1) Distribution of surface water During the last 25 years, there have been many approaches to extract water bodies from multispectral imagery (Table 1 &Table 2). These approaches can be divided into three levels: pixelwise classification, object-based image analysis (OBIA) and subpixel fraction estimation. Pixelwise approaches directly extract the pixels associated with water body areas, mainly considering the spectral characteristics of targets. OBIA groups the adjacent pixels with similar features into homogeneous clusters, which provides valuable information, including spectral, textural, shape, and spatial relationships. Subpixel fraction estimation considers the mixed pixels in remote sensing images and estimates the fraction of water bodies in each pixel using the spectral mixture analysis (SMA) method. However, the relevant methodologies fall into two categories depending on whether training samples are needed. The sample-based approach relies on the training dataset for supervised classification at both the pixel and object levels or on pure endmember selection to derive the subpixel water body fraction. The rule-based approach is based on prior knowledge of the target and background instead of the known sample data.

Small water bodies

Global water body map Urban water bodies [START_REF] Feng | A global, high-resolution (30-m) inland water body dataset for 2000: first results of a topographic-spectral classification algorithm[END_REF] Refined by constrained energy minimization Terrain metrics and prior coarse-resolution water masks

Global inland water

Object/clusterbased image analysis [START_REF] Chen | Extraction of Glacial Lake Outlines in Tibet Plateau Using Landsat 8 Imagery and Google Earth Engine[END_REF]) [START_REF] Mitkari | Extraction of Glacial Lakes in Gangotri Glacier Using Object-Based Image Analysis[END_REF]) [START_REF] Sivanpillai | Improvements in mapping water bodies using ASTER data[END_REF] (H. Xie et al., 2016) [START_REF] Zhang | Extraction of coastline in aquaculture coast from multispectral remote sensing images: Object-based region growing integrating edge detection[END_REF] Segment images using an active contour model

Combined with band indices

Unsupervised ISODATA algorithm to generate clusters k-means/ISODATA to generate clusters Object-oriented image analysis and edge detection Glacial lake outlines

Glacial lakes

Water bodies Lake and river Coastlines Among the rule-based approaches, water indices and binary thresholding-based methods are characterized by easy implementation and a high calculation efficiency and are thus widely utilized (Table 2), especially in large-scale and time series analyses [START_REF] Campos | Normalized difference water indexes have dissimilar performances in detecting seasonal and permanent water in the Sahara-Sahel transition zone[END_REF][START_REF] Pekel | High-resolution mapping of global surface water and its long-term changes[END_REF][START_REF] Thomas | Mapping inundation in the heterogeneous floodplain wetlands of the Macquarie Marshes, using Landsat Thematic Mapper[END_REF][START_REF] Zou | Continued decrease of open surface water body area in Oklahoma during 1984-2015[END_REF]. Water indices differentiate the water bodies (normally with positive values) from the background (tending to negative values). Many water indices have been designed to enhance the separation between water bodies and other land cover types (Table 2). Water indices can be directly utilized for water mapping based on thresholding [START_REF] Allen | Global extent of rivers and streams[END_REF][START_REF] Fan | Shoreline dynamics of the active Yellow River delta since the implementation of Water-Sediment Regulation Scheme: A remote-sensing and statistics-based approach[END_REF] and serve as the basis for other algorithms, such as machine learning [START_REF] Isikdogan | Surface water mapping by deep learning[END_REF], object-level segmentation [START_REF] Mitkari | Extraction of Glacial Lakes in Gangotri Glacier Using Object-Based Image Analysis[END_REF], and subpixel mapping [START_REF] Zhou | Multiscale water body extraction in urban environments from satellite images[END_REF]. [START_REF] Gao | NDWI A Normalized Difference Water Index for Remote Sensing of Vegetation Liquid Water From Space[END_REF] NIR, SWIR normalized difference water index (NDWI) [START_REF] Mcfeeters | The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features[END_REF] Green, NIR modified NDWI (MNDWI) [START_REF] Xu | Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery[END_REF] Green, SWIR automated water extraction index (AWEI) [START_REF] Feyisa | Automated Water Extraction Index: A new technique for surface water mapping using Landsat imagery[END_REF] Blue, Green, NIR, SWIRs multi-spectral water index (WuWI) [START_REF] Wang | Multi-SpectralWater Index (MuWI): A Native 10-m Multi-SpectralWater Index for accuratewater mapping on sentinel-2[END_REF] Blue, Green, NIR, SWIRs normalized difference mud index (NDMI) [START_REF] Bernstein | Quick atmospheric correction code: algorithm description and recent upgrades[END_REF] Narrow bands with wavelengths of 795 nm and 990 nm WI2015 [START_REF] Fisher | Comparing Landsat water index methods for automated water classification in eastern Australia[END_REF] Green, Red, NIR, SWIRs NDWI built-up index (NDWI-DB) (Li et al., 2016) Blue, SWIR Tasseled Cap Wetness (TCW) [START_REF] Crist | A TM Tasseled Cap equivalent transformation for reflectance factor data[END_REF] Blue, Green, Red, NIR, SWIRs normalized difference vegetation index (NDVI) (Zhu and Woodcock, 2012) Red, NIR NDWInm (H. Xie et al., 2016) Composed of a visible band and an infrared band enhanced water index (EWI) [START_REF] Wang | A simple enhanced water index (EWI) for percent surface water estimation using landsat data[END_REF] Green, Red, NIR, SWIR simple water index (SWI) [START_REF] Malahlela | Inland waterbody mapping: towards improving discrimination and extraction of inland surface water features[END_REF] Blue, SWIR LBV transformation [START_REF] Zhang | Surface Water Extraction from Landsat 8 OLI Imagery Using the LBV Transformation[END_REF] Green, Red, NIR, SWIR

(2) Dynamics of surface water Most previous large-scale and time series mapping studies of the surface water extent were conducted based on MODIS and Landsat datasets [START_REF] Aires | Comparison of visible and multi-satellite global inundation datasets at high-spatial resolution[END_REF], and the spatial resolution was generally 250 m or 30 m, respectively [START_REF] Yésou | Nine years of water resources monitoring over the middle reaches of the Yangtze River, with ENVISAT, MODIS, Beijing-1 time series, Altimetric data and field measurements[END_REF]. The dynamics of the surface water extent have been mapped at three different time scales using Landsat data, including interannual, seasonal/inner-annual, and multi/bi-temporal change analyses (Table 3). Currently, the active Landsat 7 and 8 satellites theoretically allow an 8-day repeat coverage of the land surface. However, on May 31, 2003, the scan line corrector (SLC), which compensates for the forward motion of Landsat 7, failed and led to the loss of an estimated 22 percent of the dataset. Additional cloud cover and haziness can reduce the number of valid images, which can result in a lack of sufficient data for monthly dynamic mapping. Thus, the potential use of the Sentinel-2 constellation in water body monitoring and dynamic analysis is a highly anticipated exploration.

Contributions

The main objective of this study is therefore to continuously monitor the monthly, quarterly and annual changes in surface water body dynamics in France at a 10-m spatial resolution using Sentinel-2 imagery.

The main contributions of this study include developing an automatic methodology and application of the method at the national scale.

First, an automatic rule-based superpixel (RBSP) approach is proposed. RBSP automatically analyzes the large-scale spatiotemporal variability and trends of surface water bodies. A new spectral index-based approach is proposed to delineate the water body extent. This method considers the different error sources in water body detection: (i) built-up areas using an automated water extraction index (AWEI), with noise associated with shadow and very-high albedo objects eliminated and (ii) natural areas using the normalized difference muddy index (NDMI) (for muddy and shallow water bodies) and the AWEI, with noise associated with vegetation ice and snow eliminated. Additionally, the proposed approach is implemented at the superpixel level. Superpixel segmentation groups the connected pixels with similar characteristics, and the automatic processing method is in accordance with the research objectives. RBSP is automatically run in the Google Earth Engine (GEE) platform for large-scale dynamic mapping.

Second, the Surface Water Dynamics in France (SWDF) product is generated with monthly, quarterly and annual surface water extents. To our knowledge, this product is the first national inland surface water map published at the monthly temporal scale and a 10-m spatial scale. RBSP is applied to monthly and quarterly images with medium composite pixels. Annual occurrence maps are further generated based on the pixel frequency of the monthly water maps. Our product is evaluated in both the dynamic and static states and at both national and local scales: (i) the monthly dynamic water maps are compared with satellite image based water masks at the national scale (JRC GSW monthly water history) and at local scales (over two lakes, i.e., Lake Der-Chantecoq and Lake Orient, as well as 200 random sampling points); (ii) the annual surface water extent is validated by satellite image-based water maps (JRC GSW annual product and Theia land cover product [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]) and an official database (BD Carthage hydrological product, main watercourses, and major lakes and reservoir) at the national scale and small water bodies (ponds) at local scales (including the Sologne region of the Department of Loir-et-Cher).

Methodology

The methodology part involves the overall workflow (Figure 1) of the proposed RBSP approach and the generated SWDF product. First, Sentinel-2 Level-1C data covering the France (Section 2.1) are collected.

Second, these images are preprocessed (Section 2.2) to eliminate the invalid pixels covered by clouds, to divide the bu ilt-up and natural scenes and to exclude the possible terrain shadows. Third, the proposed RBSP approach is explained, including the different spectral rules for natural and built-up scenes in Section 2.3 and superpixel technique to generate the homogeneous objects in Section 2.4. Finally, the auxiliary data and metrics used to evaluate the generated SWDF product are described in Section 2.5.

Figure 1. The workflow of the proposed RBSP approach and the validation of the obtained SWDF product

Study areas

The experiments presented in this paper focus on European France (i.e. mainland France and Corsica) (Figure 2), which encompasses an area of approximately 551,695 km 2 . France lies within the northern temperate zone and mainly includes oceanic, semi-continental, Mediterranean, and mountain climates [START_REF] Fort | Landscapes and Landforms of France (World Geomorphological Landscapes)[END_REF][START_REF] Terasmaa | A quantitative assessment of the contribution of small standing water bodies to the European waterscapes -case of Estonia and France[END_REF]. France possesses a wide variety of landscapes, from coastal plains in the north and west to the mountain ranges of the Alps in the southeast, as well as the Massif Central in the south and the Pyrenees in the southwest. Additionally, the islands of Corsica lie off the Mediterranean coast.

France has an extensive river system mainly consisting of four major rivers and their tributaries, including the Seine, Loire, Garonne and Rhône Rivers, with a combined catchment area that includes over 62% of the territory [START_REF] Kristensen | Surface water quality monitoring[END_REF]. In addition, the Rhine River and some tributaries flow through the northeastern portion of the region. Knoema Corporation (2017) reported approximately 153,000 ha of inland water occupied by major rivers, lakes and reservoirs in France. Considering the water bodies larger than 1 ha, [START_REF] Bartout | L'inventaire des plans d'eau français: Outil d'une meilleure gestion des eaux de surface[END_REF] regarded 555,000 waters bodies as corresponding to a cumulated surface area of 450,000 ha.

This study focuses on inland water detection and neglects coastline monitoring. The national and administrative boundaries are from GADM, the database of global administrative areas.

Datasets and GEE platform

This project uses the open source GEE as the research platform. The GEE provides programming and graphic interfaces for scientific applications using remote sensing data. With the powerful Google cloud storage and computational hardware technologies to accelerate remotely sensed data processing [START_REF] Gorelick | Google Earth Engine: Planetary-scale geospatial analysis for everyone[END_REF][START_REF] Trianni | Scaling up to national/regional urban extent mapping using landsat data[END_REF], GEE is advantageous for large-scale mapping and time series analysis based on multispectral images, including analyses of surface water [START_REF] Pekel | A near real-time water surface detection method based on HSV transformation of MODIS multi-Spectral time series data[END_REF], forest cover [START_REF] Hansen | High-Resolution Global Maps of 21st century forest cover change[END_REF], paddy rice planting areas [START_REF] Dong | Mapping paddy rice planting area in northeastern Asia with Landsat 8 images, phenology-based algorithm and Google Earth Engine[END_REF], and settlement areas [START_REF] Huang | Mapping major land cover dynamics in Beijing using all Landsat images in Google Earth Engine[END_REF].

The Sentinel-2 mission is composed of a two-satellite system, Sentinel-2A and Sentinel-2B, which were launched on 23 June 2015 and 07 March 2017 respectively, allowing for a high revisit frequency of approximately 5 days at the equator and 2-3 days at the mid-latitudes. The high revisit frequency is important for dynamic land cover mapping and monitoring. Sentinel-2 imagery includes 13 spectral bands (Table 4) that span from the visible (VIS) and near infrared (NIR) bands to the shortwave infrared (SWIR)

bands at different spatial resolutions on the ground ranging from 10 m to 60 m [START_REF] Drusch | Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services[END_REF]. In this study, six broad bands, including the VIS, NIR, and SWIR bands, and two narrow bands, the Red Edge 3 and 4 bands, are utilized. The SWIR and Red Edge bands with a spatial resolution of 20 m were resized to 10 m by dividing each pixel into 4 pixels with the same gray value to maintain the same spatial resolution as the VIS and NIR bands. The Sentinel-2 Level-1C dataset is the standard product for top of the atmosphere (TOA) reflectance and has been completely integrated into the GEE platform. A QA60 band is embedded in the Sentinel-2 Level-1C data, where opaque and cirrus clouds are computed based on spectral criteria. In the project, the Sentinel-2 images are filtered and collected with less than 20% cloud cover. Because the current QA60 band cannot provide accurate cloud and cloud detection results, we select a low cloud cover threshold (20%) to reduce the potential of including omission errors in cloud/cloud shadow detection.

Here, we used 8955 Sentinel-2 images (based on the 20% maximum cloud criteria) and analyzed an average of 373 Sentinel-2 observations for each month. Then, the QA60 band associated with the image was used to exclude the invalid pixels in each image. These images were acquired between March 2017 and February 2019 to cover the four seasons (starting on March 1, June 1, September 1, and December 1 for spring, summer, autumn and winter, respectively) during a two-year period. The annual map covers four quarterly datasets starting on March 1.

To complete the task of Month-Of-Year (MOY) time series analysis, images were placed into collections according to the calendar month. Image composition combines spatially overlapping images into a single image that has a medium value for each band within the monthly collection. Figure 3A displays the number of valid monthly composites of 24 months, which is the occurrence of available monthly surface water maps.

The proposed RBSP approach processes these 24 monthly time series composite data between March 2017

and February 2019 and generates the SWDF product. Over half of the land surface has been covered with over 20 MOYs, and over 90% percentage of the land surface is covered with over 17 MOYs. If utilizing quarterly composites, approximately 96% percentage of the land surface is covered with all eight seasons (Figure 3B). The detailed statistics of the valid MOYs and quarterly composites can be found in the Supplemental materials (Figure S3 and Figure S4). That is, a short period of monthly analyses may result in the occurrence of data gaps but are beneficial for observing the visible changes in hydrology in a higher temporal resolution. The auxiliary data include the JRC Global Human Settlement Layer (GHSL), HydroSHEDS, and Height

Above the Nearest Drainage (HAND) datasets, all of which are available globally and have been integrated

into the GEE platform. The GHSL contain multitemporal information layers on built-up areas, as derived from Landsat image collections (1975, 1990, 2000, and 2015) [START_REF] Pesaresi | GHS built-up grid, derived from Landsat, multitemporal[END_REF]. The latest built-up layer in 2015 was utilized to separate built-up and natural areas, and water delineation issues stemmed from different error sources [START_REF] Yang | Urban surface water body detection with suppressed built-up noise based on water indices from Sentinel-2 MSI imagery[END_REF].

Notably, terrain shadows are easily misclassified as water bodies, and digital elevation models (DEMs) are widely used to exclude terrain shadow effects. Two hydrologically relevant terrain models, HAND (Donchyts et al., 2016b;[START_REF] Silveira | Height Above the Nearest Drainage -a hydrologically relevant new terrain model[END_REF] and HydroSHEDS [START_REF] Lehner | Development and validation of a global database of lakes, reservoirs and wetlands[END_REF], were utilized to eliminate mountain shadow effects. Both models are based on high-resolution SRTM elevation data and are generally used in hydrological and remote sensing applications, such as water likelihood elimination and hill shadow correction. In this study, HydroSHEDS was used to mask mountain areas with slopes greater than 5

degrees. HAND was used to mask terrain shadow areas in flat regions with a threshold of 30 (Table 5).

Rule-based water body detection

Spectral indices highlight the pixels of objects of interest from the background, and binary thresholding can be used to delineate the object areas, which benefits time series and large-scale analyses due to the ease of implementation. Water body mapping based on water indices faces error identification, and the main noise source varies with different indices (such as the NDWI, MNDWI, and AWEI) and backgrounds (mainly including vegetation, shadows, snow and built-up objects) [START_REF] Yang | Urban surface water body detection with suppressed built-up noise based on water indices from Sentinel-2 MSI imagery[END_REF]. Compared to natural and open areas, the urban environment consists of heterogeneous human-made objects and can lead to the severe overestimation of water bodies. In this study, the AWEI, a water index for urban scenes [START_REF] Feyisa | Automated Water Extraction Index: A new technique for surface water mapping using Landsat imagery[END_REF], is adopted to distinguish water bodies from the background. However, the AWEI still faces some challenges; for example, the omission error of muddy and shallow water bodies widely can be considerable in natural areas, and commission error can occur due to the existence of building shadows and high-albedo objects in urban areas. The study area is thus divided into natural/open scenes and urban/built-up scenes to address the different types of error sources. The GHSL settlement regions designated in 2015 [START_REF] Pesaresi | GHS built-up grid, derived from Landsat, multitemporal[END_REF] are adopted to separate the natural and built-up areas, which were further processed with different rules (Table 5). Natural areas may include muddy and shallow water bodies, especially after rainfall events, and these areas could be underestimated by the AWEI. The NDMI [START_REF] Bernstein | Quick atmospheric correction code: algorithm description and recent upgrades[END_REF] is often utilized to highlight muddy and shallow water pixels, and it was originally designed as a filter to exclude those pixels and improve the accuracy of quick atmospheric correction (QUAC). Thus, a mixed water index (MWI) with large AWEI and NWI values is used to reflect the surface water extent in the natural environment. Moreover, vegetation indices can suppress shadow effects [START_REF] Yamazaki | Development of a global ~90m water body map using multi-temporal Landsat images[END_REF] and vegetation misclassification [START_REF] Zou | Continued decrease of open surface water body area in Oklahoma during 1984-2015[END_REF] in water body extraction. In this study, a mixed vegetation index (MVI) is used to eliminate errors in the surface water extent, and only those pixels that meet the criteria %&' -%(' ≤ 0.1 are classified as open surface water body pixels. Additionally, ice and snow generally display a very high degree of reflection at visible wavelengths and low reflection in the NIR and SWIR band. That is, ice and snow have a similar spectral trend (from VIS to NIR and SWIR) as water bodies, except for its stronger reflectance in the VIS bands. Thus, the blue band (greater than 0.5) is used to exclude ice and snow cover in mountainous areas.

Urban scenes involve heterogeneous human-made objects, where some building shadow areas and veryhigh albedo objects may return high positive values after AWEI calculations and serve as the main noise sources for water body maps. The NIR band is used to eliminate very-high albedo objects considering the ultralow reflectance of water bodies. Although the AWEI can suppress low-albedo objects and shadows in urban scenes, the misclassification of shadow areas still cannot be avoided. Thus, we designed an urban shadow index (USI) to highlight the shadows in urban scenes. The index considers the low reflectance of shadow areas in VIS bands compared to that of water bodies, and the reflectance values of both types of objects tends to zero in the NIR and SWIR bands.

The binary segmentation threshold is an important factor when using spectral indices to delineate target objects. The imaging of spectral indices involves polarization, where the pixel values of an object of interest tend to be positive and the background returns negative values in theory. The histogram of an index image is thus characteristic of a bimodal distribution representing the object and background and a deep and sharp valley between two peaks. Although a user-defined threshold can return optimum results, it is more appropriate to use an automated threshold or an empirical threshold to automate the delineation process, especially for large-scale and time series analyses.

In this paper, an improved Otsu threshold is adopted based on the Canny edge detection algorithm (Donchyts et al., 2016a). The performance of the global thresholding techniques (including Otsu's method) used for binary segmentation is limited for small objects and images with abundant noise [START_REF] Lee | A comparative performance study of several global thresholding techniques for segmentation[END_REF].

Land surface water represents a small fraction of the land cover in some administrative divisions in France.

Thus, the Canny edge detection algorithm is first used to identify the pixels within the buffer areas of each edge. A histogram-based Otsu approach is then applied to these pixels to filter low-probability water bodies in the scene. In addition, several soft empirical thresholds are adopted to exclude error sources to some extent. Although some pixels of these land cover types (such as snow, very-high albedo objects and vegetation) may return values similar to those of water bodies after water index calculation, they exhibit differences in other spectral bands and indices. Thus, soft empirical thresholds can effectively reduce overestimation issues. Table 5 lists the specific thresholds for different objects. For example, very-high albedo masking is given a loose threshold of 0.2 because water bodies absorb most of the spectral energy, whereas high-albedo objects have strong reflectance.

Superpixel water body mapping

Due to the lack of consideration of contextual information, the water body thematic maps obtained at the pixelwise level often experience a "salt and pepper" problem with sparse noise [START_REF] Yang | Evaluation of automated urban surface water extraction from Sentinel-2A imagery using different water indices[END_REF][START_REF] Zhang | Superpixel-based graphical model for remote sensing image mapping[END_REF]. Rather than concentrating on individual pixels, OBIA groups the nearby pixels with similar characteristics as homogeneous clusters [START_REF] Fernández | Influence of data source and training size on impervious surface areas classification using VHR satellite and aerial imagery through an objectbased approach[END_REF][START_REF] Huang | Combining Pixel-and Object-Based Machine Learning for Identification of Water-Body Types from Urban High-Resolution Remote-Sensing Imagery[END_REF][START_REF] Mitkari | Extraction of Glacial Lakes in Gangotri Glacier Using Object-Based Image Analysis[END_REF], which can effectively restrain the "salt and pepper" phenomenon and convey valuable information, including the spectral, textural, shape, and spatial information associated with adjacent objects. However, OBIA can be time consuming when extracting information from large images because most segmentation algorithms use the pixel grid as the initial object representation [START_REF] Stutz | Superpixels: An evaluation of the state-of-the-art[END_REF]. Additionally, automatic segmentation remains an unresolved problem because segmentation is sensitive to many factors, such as the image sensor resolution, scene complexity and number of bands [START_REF] Csillik | Fast segmentation and classification of very high resolution remote sensing data using SLIC superpixels[END_REF][START_REF] Gong | Superpixel-Based Difference Representation Learning for Change Detection in Multispectral Remote[END_REF].

In the field of computer science, superpixel segmentation has become increasingly popular, and images are divided into hundreds of non-overlapping superpixels [START_REF] Ren | Learning a Classification Model for Segmentation[END_REF]. Similar to OBIA, a superpixel is composed of pixels that are spectrally similar and spatially adjacent, and a superpixel is the basic unit for subsequent processing steps. Compared with OBIA, superpixel segmentation can be performed automatically for large-scale remote sensing images with low memory requirements and few parameters.

In the proposed framework, a simple non-iterative clustering (SNIC) algorithm [START_REF] Achanta | Superpixels and polygons using simple non-iterative clustering[END_REF]) is applied to generate the corresponding superpixel blocks. The SNIC algorithm is an improved version of the simple linear iterative clustering (SLIC) algorithm [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF] and is a fast and powerful algorithm with high boundary adherence and low complexity [START_REF] Gharibbafghi | Modified Superpixel Segmentation for Digital Surface Model Refinement and Building Extraction from Satellite Stereo Imagery[END_REF]. The SNIC algorithm has been proven to perform better and faster than other state-of-the-art superpixel algorithms with less memory in comparisons based on segmentation benchmarks [START_REF] Achanta | Superpixels and polygons using simple non-iterative clustering[END_REF]. This study aims to automatically and rapidly monitor water body dynamics. SNIC is thus selected due to its speed and ability to perform well for various study sites with a single set of default parameters.

The main parameter of SNIC is compactness. A large value of compactness reflects clusters with a generally rectangular shape. In this study, a low compactness value (0.1) is selected considering the irregular shape of surface water bodies. Superpixels are composed of clustered individual pixels (Figure 4), and the corresponding mean value can be used as an input for the rule-based superpixel water body delineation algorithm (Table 5). This approach results in a fast, simple, and efficient computation.

Validation and comparison

A statistically rigorous validation for this product would be desirable. However, a statistically robust national validation dataset is not available to measure the accuracy of this national water body database, especially considering the extent of surface water dynamics. Nonetheless, we evaluated our product with other existing global, national and local datasets (see S.2 in the Supplemental materials for the available online addresses). We also performed a qualitative assessment based on careful visual interpretation (see S.3 in the Supplemental materials and attached Data). Table 6 lists the main reference datasets from HR satellite image-derived water masks and official databases and the evaluation metrics.

A detection rate (+) is used to measure how well the extent (for lake) or length (for river) is accurately extracted. + is the ratio of the correctly extracted extent (for lake) or watercourse (for river) to the corresponding value in the reference dataset (Eq. 2). The extracted watercourse is calculated through an intersection process with the buffering area of the experimental river extent considering the positional deviation of the reference dataset.

+ =

.6789:; <6=</ ?9@@=6 <6=< ∩B=@=6=C:= D6=</E=CF; B=@=6=C:= <6=</E=CF;

(2)

The monthly surface water area model is a function of the sines and cosines shown in Eq. ( 3), and it is initially used to predict the day-of-year (DOY) time series surface reflectance for Landsat data (Zhu et al., 2012;[START_REF] Zhu | Continuous change detection and classification of land cover using all available Landsat data[END_REF]. In this study, we utilize the monthly composite image to simulate the month-of-year (MOY) time series of the surface area. The sine and cosine models estimate the interannual seasonal changes and inner-annual trend simultaneously, and these variations are in accord with the land surface water dynamics. A few coefficients are required to fit the functions considering the relatively small number of monthly composites. The generalized reduced gradient (GRG) solution method [START_REF] Lasdon | Design and Testing of a Generalized Reduced Gradient Code for Nonlinear Programming[END_REF] is used to fit the nonlinear time series model. Additionally, a linear regression and the correlation coefficient were used to evaluate the spatial distribution (in the units of 96 administrative divisions) of the surface water and the general trend (monthly change) of the surface water area. The correlation (a value between -1 and +1) is a common numerical measure of the degree of similarity or linear association between two variables. In this study, the Pearson correlation coefficient (\) (Eq. 4) was adopted to measure the consistency of the extent area prediction based on two products. The higher the positive value is, the more similar between the two results are in the estimation of the surface water area.

G H ∆J

\ = ∑ ^_4^̅ a _ 4a̅ bc _de f∑ ^_4^̅ g ∑ a _ 4a̅ g bc _de bc _de (4)
where h is the administrative division in France, and i j and \ j are the surface water areas within the division estimated by our SWDF product and previous products (GSW and OSO), respectively.

(2) National scale datasets Two HR satellite derived water masks and some official datasets at the national scale are available to evaluate the SWDF results. JRC GSW used 3,865,618 scenes (till now) from Landsat 5, 7, and 8 to quantify global water body dynamics from March 1984 to December 2018, at a 30-m spatial resolution, with an overall accuracy over 90%. GSW is currently the only product providing monthly dynamic data for inland surface water at a spatial resolution of several decameters. Thus, a detailed comparison of the monthly dynamics and annual surface water map in 2017 and 2018 from GSW and SWDF is conducted in this project.

The French Theia Land Data Centre has set up a Scientific Expertise Centre OSO ("Occupation des sols" in French) group, and the aim is to produce a land cover map of France using Sentinel-2 images. The OSO product is updated once a year, with the inland surface water extent at a 10-m resolution; the annual surface water extent from 2016-2018 has been released [START_REF] Inglada | Operational High Resolution Land Cover Map Production at the Country Scale Using Satellite Image Time Series[END_REF]. The overall accuracy of land classification is approximately 90%, and the F-score for the surface water is approximately 0.99. The consistency of annual water maps in 2017 and 2018 from SWDF and OSO are compared by using a confusion matrix and correlation analysis.

The national hydrological surface of the BD Carthage database provided by IGN, the French National Institute of Geography, drew the surface water extent in 2016. BD Carthage is used to evaluate the maximum annual water extents in 2017 and 2018 in the SWDF product using confusion matrix-based metrics. However, the different years between BD Carthage and SWDF mean that the surface water may vary in terms spatial distribution.

Other hydrological datasets are also used to evaluate the extraction rates (+) of the main rivers and lakes in the SWDF annual surface water maps (Figure 5). The watercourses of the major rivers in Europe, with a catchment area larger than 5,000 km 2 , released by European Environmental Agency, and the surface extents of the main lakes and reservoirs published by Système d'Information sur l'Eau (SIE) are used as the reference maps. (3) Local scale datasets Three local approaches have been carried out. The first approach corresponds to the analysis of 200 samples points' analysis, the second approach compares the monthly surface water dynamics with a detailed Sentinel-2 time series over the two major reservoirs of Lake Der-Chantecoq and Lake Orient (Figure 5), and the third approach considers a large number of ponds in the Sologne region of the Loir-et-Cher Department.

The detailed information and vector datasets of the sample points and these two reservoirs are provided in the Supplemental materials (S.3) and attached Data files.

To validate the monthly time series of surface water dynamics, we estimated confusion matrices based on 200 sample points (green points in Figure 5), which were selected from a stratified random sampling design using water/non-water strata and monthly time series of Sentinel-2 composite images. A pixel could alternate among water, non-water and no-data over the time series, and a careful visual interpretation is likely the most stable approach other than long-term field work. The accuracy of the surface water time series was described by summarizing the data in a confusion matrix and estimating the normalized coefficients (Eq. 1) of water/non-water samples across France.

The reservoirs, Lake Der-Chantecoq, 48 km 2 , and Lake Orient, 23 km 2 , (enlarged views in Figure 5) are the largest and third-largest artificial lakes in France, respectively. The Lake Der-Chantecoq and Lake Orient reservoirs are designed to protect Paris from floods by holding the water of the Marne River and Seine River, respectively. These reservoirs are fully controlled; infilling occurs from November to June, and at that time, water is taken from the Marne River (for Lake Der-Chantecoq) and the Seine (Lake Orient).

From July to October, water is released from the reservoirs for flow replenishment of the rivers. Surface water changes dramatically during the year, reducing by half of the original extent from the high to low period for Lake Der-Chantecoq and reducing to a fourth of the original extent for Lake Orient. During the period of March 2017 to February 2019, 46 and 49 Sentinel-2 images were selected over Lake Der-Chantecoq and Lake Orient, respectively. One Landsat-8 image acquired on December 7, 2017, is used because of no clear Sentinel-2 image is available. We generated the reference data of these two lakes using an SVM classification approach and an on-screen quality check.

Additionally, the Sologne region, which is a wet and relatively wild area located in the southern portion of Loir-et-Cher Department, is characterized by a large number of ponds, with a total of approximately 3,200 pounds, for a total surface area of 12,000 ha. These ponds correspond to small water bodies, where only 50 water bodies reach an area of 50 ha and the largest pond has an area of 180 ha. The reference dataset for these ponds released by the Departmental Direction of the Territories (DDT) is used to evaluate the accuracy of the proposed approach for small water body detection using the confusion matrix.

Results

Monthly and quarterly dynamics of the surface water extent

The development of a consistent and automated RBSP workflow enables us to generate a national-scale SWDF product (vector format products are attached in the Supplemental materials and Data files). We generate the spatiotemporal dynamics of inland surface water bodies using monthly composite images extending from March 2017 to February 2019. Figure 6 displays the land surface water extent in France in February 2019. The detailed results include the urban areas in Bordeaux, the natural environment of forest and farmland near Paris and in Loir-et-Cher, and the mountainous area by the Alps. Correspondingly, the land surface water bodies involve the Garonne River, which runs through Bordeaux; the Seine River, which flows through natural areas; small water bodies, such as ponds, that are surrounded by farmlands; and Lake Annecy, which is located in alpine area with terrain shadow and snow. These different types of water bodies and backgrounds reflect the robustness of the RBSP and the effectiveness of the SWDF analysis method. As the only two available products recording the monthly variation in surface water extent, SWDF and GSW are compared in terms of the aspect of area predictions. We construct two MOY time series models based on the monthly surface water area calculated in the SWDF (at the annual scale) and GSW (between March and October) products, and then, we evaluate the predicted trends based on sine and cosine fitting (Figure 7). Both products exhibit similar seasonal trends, and the correlation coefficient of the variation is 0.940. However, SWDF provides more of the larger surface water area than GSW, and this average increase of approximately 35,000 ha is related to the gain of the spatial resolution from the 30 m of Landsat image to the 10 m of Senitnel-2 image. The SWDF products describe the monthly variation based on the monthly composite data. In fact, quarterly variation can be generated if using quarterly composite data, which have also been provided in the attached Data as well. The area surface water extent in four seasons is estimated as well (Figure 8). The results indicate a general seasonality of dry autumn and wet winter. During the winter of 2017-2018 (Figure 8), particularly in January 2018 (Figure 7), France witnessed an intense episode of flooding with several successive flood waves. The highest values of SWDF MOY are related to this long flood period. To illustrate the temporal resolution of SWDF, we focused on a detailed example (Figure 9) of mapped surface water in the southern coastal areas covering parts of the Camargue Regional Nature Park, the Pond of Vaccarès, and the lower courses of the Rhône River and the Pond of Berre, which is largest salt water lake in France, within the department of Bouches-du-Rhône, at the monthly time scale. The MOY surface water area estimated from SWDF and GSW is also fitted with the sine and cosine functions (Eq. 4). The result displays dramatic seasonal behavior (Figure 9). In general, the monthly time series trend features the transition from dry summers beginning as early as June to wet winters beginning in October, which is in Here, at a more detailed scale, an example of lake dynamics across Lake Der-Chantecoq and Lake

Orient in 2017 and 2018 from the quarterly maps is shown (Figure 10). The surface water extent displays seasonal variations according to the quarterly time series analysis. Compared to the monthly variation in Figure 11, the quarterly maps are more robust because clear images are generally available for the whole of France during the quarterly time span. Nevertheless, the greater accuracy of the spatial distribution reduces the temporal variation. For example, the changes between December and February (within the period of winter) and between September and October (within the period of autumn) are noteworthy; yet, the quarterly map cannot capture these monthly changes. In the monthly product (Figure 11), Lake Der-Chantecoq is completely filled from December to February, and Lake Orient ran dry during the September and October.

The quarterly product (Figure 10) only displays the dry situation in autumn and the waterlogged state in winter. The correlation between the SWDF monthly time series products and the reference maps is 0.946 for Lake Der-Chantecoq and 0.892 for Lake Orient. That is, the monthly map can monitor the variation at higher temporal resolution with an acceptable extraction accuracy. In addition, the dramatically smallest surfaces' period occurs between November to February, which cannot be reflected in the GSW product owing to the data deficiency. Additionally, 200 sample points are marked as water bodies, non-water bodies and no data by visual interpretation based on the reference monthly composite images. In total, 4800 sample points over 24 months were obtained with two dimensions ("actual" by manual judgment and "predicted" by the RBSP method). The overall confusion matrix is presented in Table 7 after excluding 840 no-data sample points (Table S3 shows the monthly confusion matrix in the Supplemental materials). The ACC and MCC for the general products within a two-year period are 0.932 and 0.865, respectively. Owing to the wet winter season in November and December, the results during these two months present relatively low accuracy; the ACC and MCC values are approximately 0.87 and 0.75, respectively.

Water frequency and annual surface extent

The annual frequency of surface water, expressed as a value from 0 to 100%, reflects the number of times a pixel is flagged as water over the total number of cloud-free monthly compositions during the year.

Figure 12 displays the annual frequency of surface water in 2018. The satellite-based frequency of water bodies can be decreased by many factors, including the clouds (omitted by the cloud-screening algorithm), geometric mismatch and artifacts over certain areas [START_REF] Zou | Divergent trends of open-surface water body area in the contiguous United States from 1984 to 2016[END_REF]. In the project, an annual frequency of no lower than 70% is regarded as a satellite-based permanent water body. Frequencies between 20% and 70% can be regarded as seasonal water bodies, where have water some time in a year. An annual frequency of no higher than 20% indicates that the pixels are marked as water bodies during one or two months of the year.

These pixels could be the noise from cloud shadows and inundation zones, which are somewhat mixed in the annual frequency map and need further consideration. These permanent and maximum surface water maps were quantitatively compared with the HR satellite derived annual maps from 2017 and 2018 (Table 9). The confusion matrix was calculated to determine the degree of consistency of the surface water extents predicted by the different products. TP and TN denote the same extracted regions of water bodies and non-water bodies, and FP and FN denote the omission in one product and commission in another product. Additionally, correlation analysis is performed based on the surface water areas in 96 administrative divisions. The coefficient of determination (R 2 ) (Figure 13) of linear regression is found to be approximately 0.98 for GSW and 0.80 on average for OSO compared with SWDF.

The Pearson correlation coefficient (\) between the SWDF and GSW methods is approximately 0.99 on average. In general, the SWDF and GSW approaches display a high degree of consistency in both areabased correlation analysis and region-based confusion matrix aspects, perhaps due to their similar solutions based on the annual percentage of pixels.

Also, the SWDF annual products are evaluated using the BD Carthage hydrological database obtained in 2016. The confusion matrix-based metrics were 0.617 (PA), 0.801 (UA), 0.995 (ACC) and 0.700 (MCC) on average for the maximum extent and 0.611 (PA), 0.842 (UA), 0.996 (ACC) and 0.715 (MCC) for permanent water bodies. The relative high UA means that SWDF soundly overcomes the error extraction of water bodies soundly. That is, the noise is effectively excluded based on the water occurrence analysis. The low PA means that the SWDF faces greater challenges in detecting the missing parts of the surface water bodies.

However, the hydrological data were obtained before 2016 and thus the dissimilarity cannot be avoided between the different temporal data. Furthermore, the extraction rates (Figure 14) of the main rivers and main lakes are 0.929 and 0.802 in the annual water map of 2018, respectively. Further, the extraction of small water bodies is analyzed based on a large number of ponds spreading across Loir-et-Cher. The confusion matrix analysis of these small water bodies (Figure 15A) indicates that PA, UA, ACC and MCC are 0.683, 0.882, 0.995 and 0.774, respectively. These metrics for the GSW method (Figure 15B) are 0.498, 0.865, 0.993, and 0.653. Both products predict a similar degree of correctness of approximately 90% and commission errors of approximately 10%; however, they face problems with small water bodies related to low completeness and high omission errors. The SWDF method at a 10-m spatial resolution has a higher extraction rate for small water bodies and ponds in Loir-et-Cher than does the GSW method at a 30-m spatial resolution. The reference data include 2498 ponds larger than 1 ha and a total area of 8286 ha (Figure 15C). SWDF detects 1900, 2498 and 2779 ponds larger than 1 ha, 0.6 ha and 0.5 ha, respectively. GSW detects 1265, 2100 and 2498 ponds larger than 1 ha, 0.5 ha and 0.3 ha, respectively. The detailed illustration in Figure 15 indicates that the boundaries of the ponds are subject to extremely omission phenomena (redaa color), especially for GSW, which are mainly caused by the vegetation and vegetation shadows. This kind of omission error reduces the surface water area for each pond. In addition, GSW includes approximately 2000 ponds between 0.2 ha and 0.5 ha, which are generally sparse noise owing to the pixelwise mapping. The official database of the surface water extent faces the problem of the temporal changes, which would reduce its reliability in measuring SWDF dynamic products.

Discussion and perspectives

Monthly time series monitoring of surface water bodies at a 10-m resolution was performed for France between March 2017 and February 2019. The SWDF results indicate that Sentinel-2 data can provide higher temporal and spatial resolution information compared with the existing surface water extent products. The proposed RBSP approach tested here allows the frequent updating of the product based on a newly available Sentinel-2 monthly composite image. However, the proposed approach can be further improved to obtain more accurate and complete products, including the commission error of moist soil, the omission error of streams, creeks and frozen water bodies, and the data deficient because of cloud cover.

Error sources and potential improvements

Clouds and cloud shadows can cause both the omission error (Figure 16A) and commission errors (Figure 16D), and they reduce the coverage of valid observations, especially during the wet winter in France.

In this study, we utilized the Sentinel-2 QA60 band to mask clouds in images with cloud cover percentages of less than 20%. The low cloud cover threshold does not affect the proposed RBSP approach. However, the coverage of the SWDF product could be further improved if an enhanced cloud and cloud shadow algorithm are available and introduced in the preprocessing stage. Recently, the Fmask 4.0 algorithm [START_REF] Qiu | Fmask 4.0: Improved cloud and cloud shadow detection in Landsats 4-8 and Sentinel-2 imagery[END_REF], which robustly detected clouds and cloud shadows in Landsat data, was introduced into Sentinel-2. This enhanced masking algorithm is expected to (i) filter the clear pixels to restrain the commission error and (ii) add experimental data with a loose cloud cover flag to obtain a higher coverage of monthly surface water maps.

Wet soil regions distributed in farmland (through irrigation) and wetland areas are somewhat misclassified as surface water bodies based on the rules of spectral indices (Figure 16E). The division of surface water, dry land and wet land is important for monitoring the inner-annual changes in intermittent rivers and ephemeral streams. Thus, some potential ideas involve the synergistically use of Sentinel-1 SAR data [START_REF] Bousbih | Soil moisture and irrigation mapping in a semi-arid region, based on the synergetic use of Sentinel-1 and Sentinel-2 data[END_REF] and Landsat 8 thermal information [START_REF] Sadeghi | The optical trapezoid model: A novel approach to remote sensing of soil moisture applied to Sentinel-2 and Landsat-8 observations[END_REF].

RBSP typically faces omission error issues along narrow rivers (Figure 16B) because of 10 meters spatial resolution and the use of superpixel technique, especially when areas are sheltered by vegetation. It is important to restore the completeness of the watercourse in these cases. Potential development steps could be considered to restore rivers from fracture effects by (i) using prior knowledge and GIS data to obtain map layers over time and limit projection distortion, (ii) implementing salient object detection approaches and perceptual organization techniques to address the computational efficiency, and (iii) performing spectral mixture analysis to extract the endmembers of different environments.

RBSP regarded ice as the other type of land cover and excluded it from surface water maps in the current project. That is, RBSP would underestimate surface water in winter because the water bodies are frozen (Figure 16C). For the SWDF product, ice is mostly limited to mountainous lakes and hydropower reservoirs in the Alps and Pyrenees areas, which represent a very small percentage of water bodies. The monitoring of inland river and lake ice indicates significant environment and climate changes [START_REF] Yang | The past and future of global river ice[END_REF]. The detection of a surface water body can provide the basic layer for further inland ice detection and estimate the percentages and frequencies of the river and lake ice. In brief, the SWDF method may overestimate the number of water bodies in locations that contain moist soil and shaded pixels and underestimate the extent of water bodies such as streams, creeks and frozen water bodies. Commission error exists in the monthly time series maps but can be effectively excluded from the annual water maps based on frequency calculations. Such accelerated error reduces the accuracy of the SWDF flood mapping results because bi-temporal maps are applied before and after flooding. Fortunately, these errors are mainly specific to the monthly maps. The quarterly and annual surface water maps appear to be more robust than the monthly maps based on the high quality of quarterly composite images and annual frequency calculations.

Data accessibility and homogeneity of S2 and Landsat

To our knowledge, until now, JRC GSW is the only available product that provides the monthly dynamics of inland surface water bodies. However, GSW could be unavailable in certain months at certain locations. For France, GSW provides the monthly dynamics between March and October, as shown in Figure 7.

Figure 17 compares the percentage of areas with valid observations in France using Sentinel-2, Landsat 8 and GSW data (synthetic use of Landsat 7 and 8). With a high revisit frequency, Sentinel-2 (85% on average) can provide more available valid data than Landsat 8 (49% on average) and even GSW (60% on average). When limited to March through October in the dry season, the available percentage of the area is similar between the obtained product and the JRC dataset (approximately 90% on average), and the valid area decreases to 58% on average if only Landsat 8 is utilized. For the quarterly maps, Sentinel-2 provides 32 99.5% coverage, on average, in France, and Landsat 8 covers 85.6% of the entire national area on average.

The SWDF method can provide the complete seasonal dynamic changes in the inland surface water extent, although the MOY time series analysis requires further improvement to obtain complete coverage, perhaps

Conclusion

We develop an automated inland surface water detection approach and release a national surface water dynamics product. To the best of our knowledge, this study is the first to monitor the monthly dynamics of the surface water extent at a 10-m resolution over a large-scale using Sentinel-2 imagery. The RBSP algorithm runs automatically to delineate surface water bodies in different environments (such as urban scenes, agricultural fields, and mountainous areas). The SWDF product determines the water occurrence in France at monthly time steps and at a 10-m spatial resolution. The geography of France includes coastline areas, mountainous areas, plains, islands and metropolises. Thus, the RBSP approach is intuitively implemented in GEE and has the potential to generate surface water dynamics at other national scopes and even the global scope.

The obtained SWDF product is evaluated and validated based on both the HR satellite image derived water maps and official datasets. The annual surface water maps of SWDF show consistency with the publicly released land surface water maps of the JRC GSW, OSO and BD Cartage based on both the spatial distribution and surface area. These products display a high correlation coefficient of over 0.950 for surface water area prediction and high overall accuracies of over 0.995 (ACC) and approximately 0.750 (MCC) based on confusion matrix analysis. Moreover, the MOY SWDF and GSW products exhibit a similar seasonal trend, with a correlation coefficient of 0.940. The SWDF results cover all 24 months in the twoyear period, but the GSW method excludes results for the winter seasons from November to February.

Additionally, the SWDF method has a higher extraction rate than the GSW method for small water bodies due to its higher spatial resolution of 10 m. For example, the completeness of pond extractions in Loir-et-Cher is 0.683 (SWDF) and 0.498 (GSW) using the two methods, respectively. Moreover, the randomly sampled points show that the monthly water dynamics of SWDF have overall accuracies of 0.932 (ACC) and 0.865 (MCC). The detailed monitoring analysis of Lake Der-Chantecoq and Lake Orient indicates the superiority of MOY compared with conventional seasonal analysis on subtle variation monitoring.

Our ongoing work is dedicated to improving the RBSP approach and the general application of the SWDF product. The automatic RBSP approach overestimates the SWDF in wet soil and shadow areas and underestimates the areas of water bodies such as rivers, streams and frozen water bodies. An advanced cloud masking algorithm and soil moisture estimation method could be implemented to reduce the commission error. A salient edge detection approach will be considered to address the omission of river fractures. The monitoring of percentage and frequency of river and lake ice is an interesting topic to be explored in our ongoing work.

Overall, the SWDF results provide a unique opportunity for the monthly continuous mapping of the surface water extent at a 10-m scale using time series of composite images. Such monthly continuous time series of surface water dynamics benefit in-depth research on the inner-annual spatiotemporal variability in surface water changes, such as for ephemeral stream and lake monitoring, seasonal variation assessment and inundation mapping. Currently, SWDF has a valid coverage of 85% for the monthly maps and 99.5% for the quarterly maps on average. The valid coverage can be further improved if an enhance cloud and cloud shadow detection algorithm is introduced and a high cloud cover flag is used. Additionally, the homogeneous use of Sentinel-2 and Landsat 8 could further increase the valid coverage of the SWDF product.

(

  [START_REF] Du | Water bodies' mapping from Sentinel-2 imagery with Modified Normalized Difference Water Index at 10-m spatial resolution produced by sharpening the swir band[END_REF])[START_REF] Fisher | Comparing Landsat water index methods for automated water classification in eastern Australia[END_REF])(Ogilvie et al., 2018)[START_REF] Yamazaki | Development of a global ~90m water body map using multi-temporal Landsat images[END_REF] ([START_REF] Yang | Urban surface water body detection with suppressed built-up noise based on water indices from Sentinel-2 MSI imagery[END_REF]

  (i) Most works utilized Landsat series data to continuously monitor the general trend of the surface water extent over several decades.[START_REF] Rokni | Water feature extraction and change detection using multitemporal landsat imagery[END_REF] modeled the spatiotemporal changes in Lake Urmia from 2000 to 2013 and reported a dramatic decreasing trend.[START_REF] Thomas | Mapping inundation in the heterogeneous floodplain wetlands of the Macquarie Marshes, using Landsat Thematic Mapper[END_REF] mapped the inundation and flooding patterns of the Macquarie Marshes by selecting Landsat images(1989- 2010) over a range of flood magnitudes.[START_REF] Deng | Spatio-temporal change of lake water extent in Wuhan urban agglomeration based on Landsat images from 1987 to[END_REF] monitored the extent of the spatiotemporal changes of a lake in Wuhan, China, based on Landsat images from 1987 to 2015. (ii) The seasonality of surface water based on interannual monitoring has also been widely explored. The most famous work was conducted by the European Commission's Joint Research Centre (JRC)[START_REF] Pekel | High-resolution mapping of global surface water and its long-term changes[END_REF], who developed high-resolution maps of the global surface water (GSW) occurrence, change, seasonality, and transition using Landsat data at 30 meter resolution. In addition,[START_REF] Campos | Normalized difference water indexes have dissimilar performances in detecting seasonal and permanent water in the Sahara-Sahel transition zone[END_REF] derived seasonal and permanent water data between 2007 and 2011 and monitored the decrease in water resources in Africa.[START_REF] Tulbure | Spatiotemporal dynamic of surface water bodies using Landsat time-series data from 1999 to 2011[END_REF] studied spatially and temporally explicit time series of a surface water body on the Swan Coastal Plain from 1999 to 2011 and adopted landscape metrics to determine the extent of changes in the seasonally continuous surface water body after comparing summer and winter images.[START_REF] Zou | Continued decrease of open surface water body area in Oklahoma during 1984-2015[END_REF] generated four water body extent maps (maximum, year-long, seasonal, and average maps) of Oklahoma from 1984 to 2015 based on the annual water body frequency. (iii) Bi-/multitemporal images are typically applied to investigate flooding mapping and land use/cover (LULC) changes.[START_REF] Chignell | Multi-temporal independent component analysis and landsat 8 for delineating maximum extent of the 2013 Colorado front range flood[END_REF] utilized pre-and post-flood Landsat 8 images to produce a flood layer image at the regional scale of the Colorado Front Range Flood in 2013.[START_REF] Bayram | An Integrated Approach to Temporal Monitoring of the Shoreline and Basin of Terkos Lake[END_REF] analyzed the combined shoreline and LULC changes of the Terkos Lake basin using Landsat satellite imagesfrom 1986, 2001, and 2009.

Figure 2 .

 2 Figure 2. (A) The study area of France consists of 96 administrative divisions and (B) an example of a Sentinel-2 annual composite image

Figure 3 .

 3 Figure 3. Number of valid MOYs (A) and quarterly composites (B) during a two-year period using Sentinel-2 with less than 20% cloud cover

Figure 4 .

 4 Figure 4. SNIC segmentation to generate a homogeneous superpixel image using the mean value of the pixels within a block

( 1 )

 1 Evaluation metricsThe monthly and annual surface water extents are evaluated considering classification of water bodies and estimation of area extent. The classification accuracy of water body pixels is measured by the confusion matrix and detection rate (+). Except for the distribution of the water bodies, it is important to estimate the amount of the surface water bodies. Based on sine and cosine fitting, a MOY model is used to predict the tendency of the monthly surface water area. Linear regression and the correlation coefficient are also used to quantitatively analyze the consistency of the surface water area estimation.The confusion matrix divides the pixels in the study area into four classes: TP (true positive), FN (false negative), FP (false positive), and TN (true negative), reflecting accurate pixel extraction, missing water bodies, inaccurate extraction, and the accurate rejection of non-water, respectively. Four normalized metrics (Eq. 1) were then calculated to assess the performance of the proposed approach. The producer's accuracy (PA) and user's accuracy (UA) were used to indicate completeness and correctness, respectively. A low PA reflects serious omission error, and a low UA indicates an extreme commission error. The accuracy (ACC) and Matthews correlation coefficient (MCC) indicate the general accuracy of the approach. In this study, the extents of surface water and non-water bodies may be unbalanced and vary greatly. The MCC considers the four confusion matrix categories and is thus more informative than the ACC, especially when the water bodies account for a small portion of the environment.

  where ∆J is the month number of the sequence, N is the number of years of utilized Sentinel-2 data, M N is the coefficient for overall values, M O and W O are the coefficients of inner-annual change, and M T and W T are the coefficients of interannual change.

Figure 5 .

 5 Figure 5. Reference data used to evaluate the water body maps, which involve the publicly released datasets (main rivers and lakes, and ponds in Loir-et-Cher) and a detailed monitoring analysis (Lake Der-Chantecoq and Lake Orient, as well as 200 sample points)

Figure 6 .

 6 Figure 6. Surface water extent using monthly composite data during February 2019 (A) in France and several details over different environment/landscape units, including (B) urban, (C) agricultural, (D) mountainous and (E) forest areas.

Figure 7 .

 7 Figure 7. Tendency of national MOY surface water area variations estimated from the SWDF and GSW products. The crossing points indicate the monthly surface water area of France, and the curves present the fitted seasonal trend by using harmonic model. GSW provides the surface water maps between March and October.

Figure 8 .

 8 Figure 8. Surface water area of France in four seasons estimated in SWDF product

  accord with the Mediterranean climate characteristics. In addition, there are only two months (May 2017 and Nov. 2018) that are extremely affected by clouds in SWDF over the 24 MOYs. The detailed monthly changes and corresponding false color composite images are provided in Figure S1and Figure S2 (in the Supplemental materials S.4).

Figure 9 .

 9 Figure 9. Monthly time series changes in a part of the southern coastal area. (A) Monthly time series estimation of the surface water area in the study area SWDF, GSW and fitted models. (B) & (C) are two monthly maps provided by SWDF. (D-G) display two zoom-in regions, where the surface water extents (yellow lines) overlay the Sentinel-2 monthly composite image.

Figure 10 .

 10 Figure 10. Quarterly time series changes in the surface water extent of Lake Der-Chantecoq and Lake Orient. The yellow lines display the surface water extents provided by SWDF. The background corresponds to the quarterly median composite data based on Sentinel-2 images.

Figure 11 .

 11 Figure 11. Monthly time series changes in Lake Der-Chantecoq and Lake Orient. The curves (C) display the monthly surface water area estimated by SWDF and the reference values. A, B, D and E are surface water extents (yellow lines) overlaying on Sentinel-2 images. These monthly variations, including A and B in winter and D and E in autumn, cannot be reflected by quarterly maps.

Figure 12 .

 12 Figure 12. Annual surface water occurrence map in 2018 generated by calculating the frequency of pixels marked as water bodies in the monthly maps. (A) France, (B) Lake Orient, and (C) Lake Der-Chantecoq. Note: The gap inside the Lake Der-Chantecoq is from the separate processing of the two divisions (Marne and Haute-Marne).In the project, we generated water body maps in 2017 and 2018 based on the frequency map. The minimum water extent covering the permanent water bodies was 308,374 ha and 313,808 ha in 2017 and 2018, respectively, while the maximum water extent, including the permanent and seasonal water bodies,

Figure 14 .

 14 Figure 14. Evaluation of the extraction rate of main riverway and lake extents in 2018. (B) (C) and (D) are some zoomin details in different geological regions of France, which are located in Paris Basin, mountainous Alps and Mediterranean coastal area, respectively.

Figure 15 .

 15 Figure 15. Evaluation maps of pond extraction in Loir-et-Cher based on the (A) SWDF and (B) GSW products. (C) displays the statistics of the number of ponds in different size.

Figure 16 .

 16 Figure 16. Main error sources in the current SWDF products. Omission errors from (A) cloud cover, (B) stream of Claise River and (C) frozen water bodies in Plagnes Lake and Montriond Lake. Commission errors from (D) cloud shadow (E) wet land and forests shadow areas.

  by homogeneously merging Landsat and Sentinel-2 data.

Figure 3 and

 3 Figure 3 and Figure 18 present the valid MOY maps in a two-year period using Sentinel-2, Landsat 8 and Sentinel-2 and Landsat 8 combined. Our future work will involve the homogeneous merging of the Landsat and Sentinel-2 datasets, especially considering the recently published Landsat Analysis-ready Data (ARD) products and upcoming Landsat 9 scheduled for launch in December 2020.

Figure 17 .Figure 18 .

 1718 Figure 17. The percentage of the monthly composite data coverage from Sentinel-2, Landsat 8 and the combined use of Landsat 7 and Landsat 8

Table 1 .

 1 Methodologies of water body extraction from multispectral imagery

Table 2 .

 2 Water indices designed for water body detection

	Water indices	Literature	Bands
	normalized difference water index (NDWI)

Table 3 .

 3 Temporal scales for water body dynamics monitoring To our knowledge, monthly dynamics have only been assessed by[START_REF] Hui | International Journal of Remote Modelling spatialtemporal change of Poyang Lake using multitemporal Landsat imagery[END_REF],[START_REF] Dronova | Object-based analysis and change detection of major wetland cover types and their classification uncertainty during the low water period at Poyang Lake, China[END_REF],[START_REF] Campos | Normalized difference water indexes have dissimilar performances in detecting seasonal and permanent water in the Sahara-Sahel transition zone[END_REF] and[START_REF] Pekel | High-resolution mapping of global surface water and its long-term changes[END_REF].[START_REF] Hui | International Journal of Remote Modelling spatialtemporal change of Poyang Lake using multitemporal Landsat imagery[END_REF] monitored the monthly changes in Poyang Lake, China, between November 1999 and October 2000, using eight images acquired from Landsat ETM+ 5 and Landsat TM 7. Dronova et al. (2011) selected four images acquired in November 2007 and January, February, and March 2008 using the Beijing-1 microsatellite and examined the surface cover composition and its changes at Poyang Lake. Campos et al. (2012) chose one image per month between 2001 and 2010 to analyze the historical trend of the water bodies, and yet, these researchers could not monitor the real-time changes in the water bodies. JRC GSW recorded the entire history of water detection on a month-by-month basis between March 1984 and December 2018. However, the coverage of the monthly map is still limited by an insufficient number of valid observations based on Landsat data.

	Temporal scale	Literature	Methodology
	Annual to decade trend	(Allen and Pavelsky, 2018)	Water indices
		(Arvor et al., 2018)	Time series indices
		(Avisse et al., 2017)	Fmask, water and vegetation indices
		(Carroll and Loboda, 2017)	Using the DSWE product
		(Deng et al., 2017)	Indices and random forest
		(Fan et al., 2018)	Water index
		(Ogilvie et al. 2018)	Water index and hydrological modeling
		(Pardo-Pascual et al., 2012)	Shoreline subpixel detection
		(Sagar et al., 2017)	Median pixel compositing of NDWI stacks
		(Shi et al., 2017)	An 'eight-field' morphological method
		(Tseng et al., 2016)	Surface water area and level changes
		(Liu and Yue, 2017)	Band value and Otsu threshold
		(Zou et al., 2018)	The relationship between water and vegetation indices
		(Yamazaki et al., 2015)	Spectral indices
	Seasonality	(Pekel et al., 2016)	Expert systems, visual analytics, and evidential reasoning
		(Tulbure and Broich, 2013)	Decision tree classification algorithm
		(Zou et al., 2017)	Spectral indices
		(Tulbure et al., 2016)	Random forest
		(Sheng et al., 2016)	Water index
	Monthly changes	(Campos et al., 2012)

[START_REF] Hui | International Journal of Remote Modelling spatialtemporal change of Poyang Lake using multitemporal Landsat imagery[END_REF]

)

[START_REF] Pekel | High-resolution mapping of global surface water and its long-term changes[END_REF]

Water indices based on one image per month between 2001 and 2010 Water indices and empirical threshold 6

Table 4 .

 4 Band information of the Sentinel-2 Level-1C data

	Band number	Band name	Resolution (m)	Band number	Band name	Resolution (m)
	B1	Aerosols	60	B8	NIR	10
	B2	Blue	10	B8A	Red Edge 4	20

Table 5 .

 5 The image characteristics used to delineate water bodies from the background using Sentinel-2

	Scenes	Feature	Equation or Methodology	Thresholding	Objectives
		HydroSHEDS	HydroSHEDS.slope		Empirical value: 5	Terrain shadows
	Preprocessing	HAND	Combined with JRC water occurrence	Empirical value: 30	Mountain shadows in flat areas
		AWEIsh	B2 + 2.5 × B3 -1.5 × B8 + B11 -0.25 × B12	Edge-based Otsu	Water maps
	Built-up areas	NIR band	B8		Empirical threshold: 0.2	Very-high albedo noise
		USI	-2 1 -B2 × 1 -B3 × 1 -B4	Edge-based Otsu	Urban shadow noise
		MWI	max NDMI, AWEI ; NDMI =	B7 -B8A B7 + B8A	Edge-based Otsu	Water maps
	Natural areas	MVI	MVI = EVI = B8 + 6 × B4 -7.5 × B2 + 1 EVI + NDVI 2 ; NDVI = B8 -B4 B8 + B4 ; 2.5 * B8 -B4	Relative threshold %&' -%(' > 0.1	Vegetation noise
		Blue band	B2		Empirical threshold: 0.5	Ice and snow noise

Table 6 .

 6 Datasets used to evaluate SWDF monthly and annual surface water maps at national and local scales

			Monthly dynamic map		Annual static map
			Dataset	Evaluation	Dataset	Evaluation
		Satellite image based product	JRC GSW monthly water history	Correlation analysis & trend analysis	Theia OSO annual product JRC GSW annual product	Confusion matrix & correlation analysis
	National	Official database	Not available		BD Carthage in 2016 Main watercourse Main lake and reservoir	Confusion matrix Extraction rate Extraction rate
	Local	Satellite image based product	200 random sample points Two seasonal lakes South coastal area	Confusion matrix Trend analyses Qualitative analyses	Small water bodies in Loir-et-Cher provided from JRC GSW	Comparison based on confusion matrix & number count
		Official database	Not available		Small water bodies in Loir-et-Cher	Confusion matrix & Number count

Table 7

 7

		. Accuracy assessment of sample points
	840 points: no data available		Actual reference by visual interpretation
	3960 points: valid MOYs		Water	Non-water
	Predicted SWDF	Water	1589	61
	by the RBSP method	Non-water	207	2103
		PA = 0.885; UA = 0.963; ACC = 0.932; MCC = 0.865

  Table 8 listed the surface water area estimated from different HR derived products and the hydrological BD Cartage dataset.

Table 8 .

 8 Surface water area in France estimated from different products

	Surface water area (ha)	SWDF permanent maximum	permanent	GSW	maximum	OSO maximum	BD Cartage in 2016 permanent maximum
	2017 2018	308,374 313,808	421,857 414,449	277,959 287,551		394,297 402,313	541,912 540,353	428,401	543,319

Table 9 .

 9 Comparison of the annual surface water maximum and permanent extent provided by three satellite derived Correlation analysis of the annual surface water maximum extent between the SWDF and GSW/OSO methods based on the surface water areas in 96 administrative divisions in (A) 2017 and (B) 2018

			products and a hydrological database
	(2017/2018)	SWDF vs OSO	SWDF vs GSW	OSO vs GSW	SWDF vs BD Cartage
		r	0.952 / 0.838	0.991 / 0.989	0.946 / 0.856	0.915 / 0.914
	Maximum extent	ACC	0.996 / 0.995	0.997 / 0.997	0.995 / 0.994	0.995 / 0.995
		MCC	0.750 / 0.700	0.779 / 0.776	0.740 / 0.687	0.698 / 0.702
	Permanent	ACC	Not available	0.997 / 0.997	Not available	0.996 / 0.996
		MCC	Not available	0.757 / 0.759	Not available	0.721 / 0.710

Acknowledgments

This work is funded by National Key R&D Program Funding (2017YFB0503905-05). The work of X.

Yang was supported by the China Scholarship Council (No. 201504490008). We thank the European Commission's Joint Research Centre, Theia, the European Environment Agency, and Système d'Information sur l'Eau (SIE) for providing the reference maps used to validate our products. The Lake de Der and Lake d Orient monitoring was realized within the framework of the SWOT Project.