Automatic Inundation Mapping Using Sentinel-2 Data Applicable to Both Camargue and Doñana Biosphere Reserves

: Flooding periodicity is crucial for biomass production and ecosystem functions in wetland areas. Local monitoring networks may be enriched by spaceborne derived products with a temporal resolution of a few days. Unsupervised computer vision techniques are preferred, since human interference and the use of training data may be kept to a minimum. Recently, a novel automatic local thresholding unsupervised methodology for separating inundated areas from non-inundated ones led to successful results for the Doñana Biosphere Reserve. This study examines the applicability of this approach to Camarque Biosphere Reserve, and proposes alternatives to the original approach to enhance accuracy and applicability for both Camargue and Doñana wetlands in a scientiﬁc quest for methods that may serve accurately biomes at both protected areas. In particular, it examines alternative inputs for automatically estimating thresholds while applying various algorithms for estimating the splitting thresholds. Reference maps for Camargue are provided by local authorities, and generated using Sentinel-2 Band 8A (NIR) and Band 12 (SWIR-2). The alternative approaches examined led to high inundation mapping accuracy. In particular, for the Camargue study area and 39 di ﬀ erent dates, the alternative approach with the highest overall Kappa coe ﬃ cient is 0.84, while, for the Doñana Biosphere Reserve and Doñana marshland (a subset of Doñana Reserve) and 7 di ﬀ erent dates, is 0.85 and 0.94, respectively. Moreover, there are alternative approaches with high overall Kappa for all areas, i.e., at 0.79 for Camargue, over 0.91 for Doñana marshland, and over 0.82 for Doñana Reserve. Additionally, this study identiﬁes the alternative approaches that perform better when the study area is extensively covered by temporary ﬂooded and emergent vegetation areas (i.e., Camargue Reserve and Doñana marshland) or when it contains a large percentage of dry areas (i.e., Doñana Reserve). The development of credible automatic thresholding techniques that can be applied to di ﬀ erent wetlands could lead to a higher degree of automation for map production, while enhancing service utilization by non-trained personnel.

be identified as possibly applicable to further wetland areas and biomes. Each examined alternative approach relies on a specific band or band combination, acknowledged as effective by the underlying physics, and a specific approach for estimating splitting thresholds. The different Sentinel-2 (S2) based inputs examined for estimating thresholds include: (i) Band 11 (SWIR-1), (ii) product (result of the multiplication) of Band 12 (SWIR-2) and Band 8A (NIR), and (iii) product of SWIR-1 and NIR. The different methods for estimating splitting thresholds include: (i) minimum entropy thresholding, and (ii) Otsu's algorithm. The results of the alternative approaches are compared against reference maps, provided for Doñana and Camargue by local research institutes, based on locally developed water detection models [22,24].

Camargue
The Camargue ( Figure 1) is a polderized delta created by the Rhône River that covers an area of 145,300 ha rarely exceeding 5 m in altitude. It comprises a high diversity of wetland types according to a water-salinity gradient such as lagoons, salt meadows, dense halophilous scrubs and steppes, brackish/freshwater marshes with tall emergent or permanent aquatic vegetation, and temporary pools. The climate is Mediterranean, being characterized by mild and wet winters and hot and dry summers. With mean annual precipitations of 600 mm, mainly concentrated from autumn to early spring, and a mean evapotranspiration of 1400 mm, the Camargue is characterized by high water deficits, especially in the summer [44,45]. As a result, 730 millions of cubic meters of water are pumped from the Rhône on average each year to compensate for river embankment, avoid soil salinization, and enhance primary production. This water, primarily pumped from March to September, is distributed through a complex network of channels for irrigating crops and pastures, as well as for flooding marshes, which support nature conservation, wildfowl hunting, and reed harvest. Accordingly, some Camargue wetlands are flooded year round or during most of the year as a result of artificial irrigation (e.g., hunting reed marsh). Other wetlands become naturally dry during the summer season (e.g., harvested reed marsh, bulrush marsh) and, lastly, a few are flooded only during a period of high rainfalls (e.g., halophilous scrubs, temporary pools) [46]. Within the context of climate change and increasing human pressures, the development of a remotely-sensed tool to monitor water seasonal variation in these wetlands is essential [47].

Doñana
The Doñana wetlands (Figure 2), covering 108,429 ha, lie within and around the delta of the Guadalquivir River in Southwest Spain and contain two main habitat types: seasonal marshes and adjacent eolian sands holding aquifer-fed dune ponds. These are surrounded by scrublands, pine forests, and cultivated areas. The marshland area is comprised of seasonal open marsh with emergent plants, temporary pools with annual plant species, and scattered halophilous scrubs [48]. The Doñana climate is Mediterranean sub-humid with hot and dry summers and mild and wet winters. The annual precipitation, with an average of 550 mm, occurs mainly between October and April and is almost absent between May and September. The highest monthly rainfall usually occurs in November and maximum water levels are reached during February. Marshes dry up slowly in late spring and most of their surface gets completely dry by the end of July [49]. Marshland's depth, turbidity, and vegetation cover varies depending on the amount and seasonal pattern of precipitation [24]. Marshes are breeding ground and host many species of migratory birds during the winter [49]. In parallel, rice-paddies, aquaculture ponds, and salt pans constitute an extensive system of artificial wetlands in Doñana that also provide a habitat for water birds and other species, especially when natural wetlands dry up seasonally [50].

Doñana
The Doñana wetlands (Figure 2), covering 108,429 ha, lie within and around the delta of the Guadalquivir River in Southwest Spain and contain two main habitat types: seasonal marshes and adjacent eolian sands holding aquifer-fed dune ponds. These are surrounded by scrublands, pine forests, and cultivated areas. The marshland area is comprised of seasonal open marsh with emergent plants, temporary pools with annual plant species, and scattered halophilous scrubs [48]. The Doñana climate is Mediterranean sub-humid with hot and dry summers and mild and wet winters. The annual precipitation, with an average of 550 mm, occurs mainly between October and April and is almost absent between May and September. The highest monthly rainfall usually occurs in November and maximum water levels are reached during February. Marshes dry up slowly in late spring and most of their surface gets completely dry by the end of July [49]. Marshland's depth, turbidity, and vegetation cover varies depending on the amount and seasonal pattern of precipitation [24]. Marshes are breeding ground and host many species of migratory birds during the winter [49]. In parallel, rice-paddies, aquaculture ponds, and salt pans constitute an extensive system of artificial wetlands in Doñana that also provide a habitat for water birds and other species, especially when natural wetlands dry up seasonally [50].

Camargue
Thirty-nine cloud-free and atmospherically corrected S2 Level-2A (L2A) products of the Camargue from 12 June, 2017 to 19 June, 2018 were downloaded from the Copernicus European Space Agency (ESA) hub (see Table 1). The shapefile containing the boundaries of Camargue comprises two different tiles (namely 31TEJ, 31TFJ). The tile products are mosaicked and clipped to the extent of the shapefile. The reference maps for Camargue were obtained by dichotomous partitioning of reflectance values encoded as 1 for water presence and 0 for water absence based on ground-truth (n = 1229) and optical-space derived (n = 2603) reference points covering the whole Biosphere Reserve area and all the main habitat types [22]. Ground-truth data refer to water level measures in different wetland types, focusing on those with a dense vegetation cover. For optical-space data, a water formula developed with SPOT-5 (Satellite Pour l'Observation de la Terre -5) [13] was used to generate a water mask for dates when Landsat-8 scenes were also available. Fifty random points from each of the main land cover types in the Camargue (n = 17) were selected on the SPOT-5 water maps and transferred to Landsat-8 scenes with similar dates. These points were then used as training data along with ground-truth measures for creating a Landsat-8 water formula. The same procedure was repeated to develop a water formula with S2 from Landsat-8. Data mining was performed with the Rpart package in the R software using eight bands of S2, as well as various spectral indices found in the literature for explaining the variables (for more details, see References [13,22]). A random selection of 30% of all points was excluded from the sample and used for validation. The best model selected used the near (NIR) and short-wave (SWIR-2) infrared wavelengths. NIR was useful for discriminating areas that are completely dry, while SWIR-2 was efficient for water detection. This model provided an overall accuracy of 94% for predicting water presence/absence with Kappa coefficients of 0.82 on both the training and validation samples [22]. Wetlands characterized by a dense cover of vegetation were correctly classified at 89%.

Doñana
Seven cloud-free and atmospherically corrected S2 L2A products of the Doñana from 1 June, 2017 to 17 April, 2018 were downloaded from the ESA hub, so that they timely overlap with Landsat cloud-free parallel data (see Table 2). The shapefile containing the boundaries of Doñana Biosphere Reserve comprises three different tiles (namely 29SQA, 29SQB, and 29SPB). The tile products were mosaicked and clipped to the extent of the shapefile. Inundation maps of the Doñana area with 30-m pixel resolution, which are generated from Landsat satellite data and provided by the Doñana Biological Station (EBD), are used as ground truth data. These reference maps were obtained by dichotomous partitioning of reflectance values on Landsat Thematic Mapper (TM) and Enhanced TM (ETM) band 5 (SWIR) based on 6005 ground-truth field sampling points, mostly collected within the marshland area of Doñana [24]. This model provided an overall accuracy of 93% for predicting water presence/absence with a Kappa coefficient of 0.65.

Methodology
The work presented in Reference [23] introduces an unsupervised approach, detecting automatically thresholds on the SWIR-1 band and on a Modified-Normalized Difference Vegetation Index (MNDVI) (estimated as the normalized difference of Band 7 and Band 5 of S2), to estimate open-water and water-vegetation subclasses. This approach demonstrated high classification accuracy for Doñana, with an overall Kappa coefficient reaching 0.94 and 0.88 for the marshland and the complete area (i.e., the Biosphere Reserve), respectively. This paper examines alternatives of the original thresholding approach driven by the findings for Camargue. Instead of the SWIR-1 band, other algebraic band combinations are also examined. Additionally, both the minimum cross entropy thresholding algorithm (MCET) [51] and Otsu's algorithm [28] are used for estimating thresholds partitioning inundated and non-inundated pixels inside an image subset based on their class distribution. The methodology's steps are presented in Figure 3.

Segmentation of the Satellite Image
As the first step, the satellite image is segmented into non-overlapping regions by utilizing the mean-shift segmentation algorithm [52]. The resulting segmentation map is utilized for selecting segments with a high percentage of inundated pixels.
The input to the segmentation algorithm is a false color image composed of Band 2 (BLUE), Band 3 (GREEN), and Band 4 (RED), which have been normalized to the range [0,255], by relying on minimum and maximum values representing the intensity percentile range from 1% to 99% per band. These bands are selected due to their 10-m resolution, which allows for a more accurate segmentation of the satellite image, compared to selection of other bands with a lower resolution [23].

Mapping of the Open-Water Subclass
Three input data alternatives are examined for estimating an initial threshold Tinit separating inundated from non-inundated pixels: (i) SWIR-1 band, which is denoted as Alt1, (ii) product (per pixel multiplication) of SWIR-2 and NIR, which is denoted as Alt2, and (iii) product of SWIR-1 band

Segmentation of the Satellite Image
As the first step, the satellite image is segmented into non-overlapping regions by utilizing the mean-shift segmentation algorithm [52]. The resulting segmentation map is utilized for selecting segments with a high percentage of inundated pixels.
The input to the segmentation algorithm is a false color image composed of Band 2 (BLUE), Band 3 (GREEN), and Band 4 (RED), which have been normalized to the range [0,255], by relying on minimum and maximum values representing the intensity percentile range from 1% to 99% per band. These bands are selected due to their 10-m resolution, which allows for a more accurate segmentation of the satellite image, compared to selection of other bands with a lower resolution [23]. Three input data alternatives are examined for estimating an initial threshold T init separating inundated from non-inundated pixels: (i) SWIR-1 band, which is denoted as Alt1, (ii) product (per pixel multiplication) of SWIR-2 and NIR, which is denoted as Alt2, and (iii) product of SWIR-1 band and NIR, which is denoted as Alt3. The first alternative Alt1 is the one proposed in the original approach [23]. The second alternative Alt2 is inspired by the approach suggested in Reference [22] for Camargue, which applies strict thresholds to SWIR-2 and NIR (see paragraph 2.2.1), and the third alternative Alt3 is a variant of the Atl2 where SWIR-1 is used instead of SWIR-2. Each of Alt1, Alt2, and Alt3 is normalized in the range of [0,255], by relying on minimum and maximum values representing the intensity percentile range from 1% to 99%. The term "Inp", used in the following, is a generic term for the input data, corresponding to either Alt1, Alt2, or Alt3.
The histogram of Inp is used to estimate an initial threshold T init separating inundated from non-inundated pixels. T init is the Inp value for which the first deep valley of this histogram is detected ( Figure 4 shows histograms using Alt2 as input). If a pixel p has Inp(p) < T init , it is denoted as inundated. Otherwise, it is denoted as non-inundated. Therefore, an initial inundation map is generated. Based on this inundation map, segments G m , where m = 1,2, . . . , M, having a large percentage of inundated pixels (i.e., over 70% of a segment's pixels are inundated) are selected and their corresponding centroids C m are estimated.
Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 20 Square patches P of expanding window size are centered around each Cm. The window size in pixels is given by [20•k × 20•k], where k = 1,2,…,20. The splitting thresholds of patches with bimodal histogram (i.e., two distinctive classes of intensity values corresponding to inundated and non-inundated classes appear in the histogram) are estimated based on the histogram thresholding algorithms (MCET or Otsu's algorithm). Their median is assumed to be the optimal threshold corresponding to the segment Gm. The usefulness of using expanding patches is to calculate the optimal threshold per segment more robustly, since its estimation is based on multiple splitting thresholds. Then, the median of selected segments' optimal thresholds is estimated as:  When MCET or Otsu is used for estimating the splitting thresholds, Tfinal is denoted as or , respectively. Average (abbreviated as "avg"): / is also examined as the final threshold. Square patches P k m of expanding window size are centered around each C m . The window size in pixels is given by [20·k × 20·k], where k = 1,2, . . . , 20. The splitting thresholds f k m of patches with bimodal histogram (i.e., two distinctive classes of intensity values corresponding to inundated and non-inundated classes appear in the histogram) are estimated based on the histogram thresholding algorithms (MCET or Otsu's algorithm). Their median is assumed to be the optimal threshold f m opt corresponding to the segment G m . The usefulness of using expanding patches is to calculate the optimal threshold per segment more robustly, since its estimation is based on multiple splitting thresholds. Then, the median of selected segments' optimal thresholds is estimated as: M opt = median Final threshold T final (see Figure 4), discriminating the open-water subclass, is estimated as max(M opt ,T init ). Pixels p with Inp(p) < T final are assumed to belong to the open-water subclass, comprising of water or water with sparse vegetation pixels.

Mapping of the Water-Vegetation Subclass
When MCET or Otsu is used for estimating the splitting thresholds, T final is denoted as T MCET final or T Otsu final , respectively. Average (abbreviated as "avg"): T avg final = T MCET final + T Otsu final /2 is also examined as the final threshold.

Mapping of the Water-Vegetation Subclass
Inundated areas may be covered by emergent vegetation. The Inp value of the pixels in these areas is higher compared to the Inp values of the pixels with open water or water with sparse vegetation belonging to the open-water class. Therefore, a threshold T upper , which is the Inp value for which the deep valley right after the first deep valley is detected in the Inp histogram, is assumed to represent the upper threshold for pixels comprising areas of water covered by dense vegetation (see Figure 4b). At the same time, areas with dense emergent vegetation are expected to have high MNDVI value. Thus, a threshold T MNDVI , which should be more than 0.4, is selected. In particular, T MNDVI is set equal to the MNDVI value, for which the first valley in the part of the histogram with MNDVI values over 0.4 is detected (see Figure 5b) following the findings in Reference [23]. When MCET or Otsu is used for estimating the splitting thresholds, Tfinal is denoted as or , respectively. Average (abbreviated as "avg"): / is also examined as the final threshold.

Comparison of Automatic Thresholding Results Against Reference Map of Camargue
The three input alternatives (Alt1 or Alt2, or Alt3) and the three different ways of estimating the final threshold (T MCET final or T Otsu final or T avg final ) form nine different alternatives, with each one combining an input alternative and a way of estimating the final threshold. The alternative ways of estimating final threshold are referred to as "MCET" (standing for T MCET final ), "OTSU" (standing for T Otsu final ), and "avg" (standing for T avg final ). For example, if "Alt2" is the input and "avg" is the way of estimating the final threshold, then the alternative will be named as "Alt2 (avg)". Regarding Camargue, each of the alternatives is applied to the 39 S2 products (see Table 1), and the results are compared to the reference maps. Pixels corresponding to the coastal waters are not taken into consideration for the comparison to obtain unbiased classification results. The results are presented using Kappa coefficient, which is estimated in terms of the observed agreement (p o ) and the agreement expected by chance (p e ) as: Kappa= (p o − p e ) / (1-p e ) [53].
The curves presented in Figure 6 show how the Kappa coefficient varies for the time period between 12 June, 2017 and 19 June, 2018 for five out of nine alternatives with the best agreement with the reference map. The best agreement is given by "Alt3 (OTSU)". "Alt3 (avg)" and "Alt2 (avg)" rank second and third, respectively, and have similar results. Lastly, "Alt3 (MCET)" and "Alt2 (MCET)" rank fourth and fifth, respectively, and give very close results. The previously mentioned ranking is based on the overall Kappa accuracy provided in Section 3.3. For all curves, it is seen that agreement decreases in the winter months from XII to II. For this period, Kappa is below 0.7 for most of the dates, while, for the rest of the year, the Kappa value is generally more than 0.7.
are met: Tfinal < Inp(p) < Tupper & MNDVI(p) > TMNDVI, presuming that Tupper and TMNDVI can be detected in the histograms of Inp and MNDVI, respectively. These latter conditions are valid for Doñana on 20 August, 2017 (see Figure 4b and Figure 5b), but not for Camargue on 18 August, 2017 (see Figure  4a and Figure 5a).

Comparison of Automatic Thresholding Results Against Reference Map of Camargue
The three input alternatives (Alt1 or Alt2, or Alt3) and the three different ways of estimating the final threshold ( or or ) form nine different alternatives, with each one combining an input alternative and a way of estimating the final threshold. The alternative ways of estimating final threshold are referred to as "MCET" (standing for ), "OTSU" (standing for ), and "avg" (standing for ). For example, if "Alt2" is the input and "avg" is the way of estimating the final threshold, then the alternative will be named as "Alt2 (avg)".
Regarding Camargue, each of the alternatives is applied to the 39 S2 products (see Table 1), and the results are compared to the reference maps. Pixels corresponding to the coastal waters are not taken into consideration for the comparison to obtain unbiased classification results. The results are presented using Kappa coefficient, which is estimated in terms of the observed agreement (po) and the agreement expected by chance (pe) as: Kappa= (po -pe) / (1-pe) [53].    (Figure 7a), compared to "Alt3 (MCET)" (Figure 7d), which largely underestimates water presence, and "Alt3 (avg)" (Figure 7c), which shows intermediate results. However, while "Atl3 (OTSU)" is able to detect water in vegetated areas more efficiently, such as the ones enclosed in the yellow squares, it tends to overestimate water presence in agricultural areas.
pink colors indicate classification differences between the reference map and the Alt3 maps. The map of "Alt3 (OTSU)" (Figure 7b) agrees better with the reference map (Figure 7a), compared to "Alt3 (MCET)" (Figure 7d), which largely underestimates water presence, and "Alt3 (avg)" (Figure 7c), which shows intermediate results. However, while "Atl3 (OTSU)" is able to detect water in vegetated areas more efficiently, such as the ones enclosed in the yellow squares, it tends to overestimate water presence in agricultural areas.

Comparison of Automatic and Thresholding Results Against Landsat Reference Maps of Doñana
In order to account for the uncertainty caused by the lower spatial resolution of the Landsat-derived reference maps, pixels in the transition zones between inundated and non-inundated areas (i.e., pixels including in their eight closest-neighbor pixels at least one pixel of a different class) in the S2 maps were excluded from the accuracy estimation (see Reference [23] for more information). Pixels in the area corresponding to sea coastal waters were also excluded. The curves, presented in Figure 8, show how Kappa coefficient varies for seven different dates regarding the five top ranking alternatives at the Doñana marshland (Figure 8a) and Doñana complete area (Figure 8b). derived reference maps, pixels in the transition zones between inundated and non-inundated areas (i.e., pixels including in their eight closest-neighbor pixels at least one pixel of a different class) in the S2 maps were excluded from the accuracy estimation (see Reference [23] for more information). Pixels in the area corresponding to sea coastal waters were also excluded. The curves, presented in Figure  8, show how Kappa coefficient varies for seven different dates regarding the five top ranking alternatives at the Doñana marshland ( Figure 8a) and Doñana complete area (Figure 8b). In Doñana marshland, "Alt1 (OTSU)" provides the highest accuracy, followed by "Alt2 (MCET)" and "Alt1 (Avg)" (see section 3.3 for more details on the ranking that is based on overall Kappa).
In the Doñana complete area, "Alt3 (MCET)" provides the highest accuracy, followed by "Alt1 (MCET)" and "Alt2 (MCET)" (see section 3.3 for more details on the ranking that is based on overall Kappa). "Alt1 (avg)" provides the most consistent Kappa values across the study period. Figure 9 shows an example of the Landsat derived inundation reference map on 27 January, 2018 and the inundation maps estimated based on the alternative approaches using Alt3 input. Green and pink colors indicate the classification differences between the reference map and the alternative maps.  In Doñana marshland, "Alt1 (OTSU)" provides the highest accuracy, followed by "Alt2 (MCET)" and "Alt1 (Avg)" (see Section 3.3 for more details on the ranking that is based on overall Kappa).
In the Doñana complete area, "Alt3 (MCET)" provides the highest accuracy, followed by "Alt1 (MCET)" and "Alt2 (MCET)" (see Section 3.3 for more details on the ranking that is based on overall Kappa). "Alt1 (avg)" provides the most consistent Kappa values across the study period. Figure 9 shows an example of the Landsat derived inundation reference map on 27 January, 2018 and the inundation maps estimated based on the alternative approaches using Alt3 input. Green and pink colors indicate the classification differences between the reference map and the alternative maps. "Alt3     Figure 10 provides overall Kappa coefficients, which are estimated when the number of true positive pixels of the inundated class, false positive pixels of the inundated class, true positive pixels of the non-inundated class and false positive pixels of the non-inundated class are added for dates per approach and examined study. The combined Overall Accuracy of each approach per study case is also given as a number close the point corresponding to its overall Kappa. Kappa values can be classified according to Reference [53] as follows: "Moderate" when 0.40 < Kappa ≤ 0.60 (shortened as "Mod"), "substantial" when 0.60 < Kappa ≤ 0.80 (shortened as "Sub") or "almost perfect" when 0.80 < Kappa ≤ 1 shortened as ("Alm Perf").

Discussion
The aim of this paper is to examine the performance of various alternative automatic thresholding approaches for inundation mapping in two different wetland areas, which include Camargue and Doñana. These are considered as testbeds to evaluate the transferability of the different approaches and select the ones that favor the acquisition of consistently credible results for both wetlands. In this study, the local thresholding alternative approaches utilize multispectral satellite data, while the vast majority of local thresholding approaches have been designed for utilizing radar data (e.g., [27,[54][55][56][57][58]). Moreover, the estimation of local thresholds relies on patches of expanding size, contrary to most of the approaches using image subsets of fixed size [55][56][57][58]. This helps to increase robustness, when input data have different spatial resolution such as when inputs from different satellites are used.
The overall Kappa of the alternative approaches varies from 0.73 to 0.84 for Camargue, from 0.74 to 0.94 for Doñana marshland, and from 0.45 to 0.85 for the complete area of Doñana. The experimental results demonstrate that there is not a common approach across all study cases achieving absolute top accuracy. "Alt3 (OTSU)" is best for Camargue, "Alt1 (OTSU)" is best for Doñana marshland, and "Alt3 (MCET)" is best for Doñana as a whole. The use of "Alt2" and "Alt3" Figure 10. Overall Kappa per approach and examined study area. The dotted line delineates the threshold of the "Alm Perf" case. Yellow arrows indicate alternative approaches succeeding "Alm Perf" results in Camargue and Doñana marshland at the same time, while blue arrows showcase the most successful overall alternative ones. The combined overall accuracy of each approach per study area is given as a number close the point corresponding to its overall Kappa. "Alt2 (avg)" and "Alt3 (avg)" (see yellow arrows in Figure 10) have both "Alm Perf" overall Kappa when considering Camargue and Doñana marshland study areas. While, for the Doñana complete area, their Kappa is "Sub", since overall Kappa for "Alt2 (avg)" and "Alt3 (avg)" is 0.63 and 0.73, respectively. On the other hand, "Alt2 (MCET)" and "Alt3 (MCET)" (see blue arrows in Figure 10) exhibit more consistent performance across sites, since their Kappa is "Alm Perf" for both Doñana marshland and Doñana complete area, while, for Camargue, their overall Kappa is very close to "Alm Perf."

Discussion
The aim of this paper is to examine the performance of various alternative automatic thresholding approaches for inundation mapping in two different wetland areas, which include Camargue and Doñana. These are considered as testbeds to evaluate the transferability of the different approaches and select the ones that favor the acquisition of consistently credible results for both wetlands. In this study, the local thresholding alternative approaches utilize multispectral satellite data, while the vast majority of local thresholding approaches have been designed for utilizing radar data (e.g., [27,[54][55][56][57][58]). Moreover, the estimation of local thresholds relies on patches of expanding size, contrary to most of the approaches using image subsets of fixed size [55][56][57][58]. This helps to increase robustness, when input data have different spatial resolution such as when inputs from different satellites are used.
The overall Kappa of the alternative approaches varies from 0.73 to 0.84 for Camargue, from 0.74 to 0.94 for Doñana marshland, and from 0.45 to 0.85 for the complete area of Doñana. The experimental results demonstrate that there is not a common approach across all study cases achieving absolute top accuracy. "Alt3 (OTSU)" is best for Camargue, "Alt1 (OTSU)" is best for Doñana marshland, and "Alt3 (MCET)" is best for Doñana as a whole. The use of "Alt2" and "Alt3" input alternatives leads to the estimation of additional inundated areas compared to "Alt1", but sometimes leads to overestimating the water presence in dry areas. "OTSU" algorithm identifies the final threshold at a relatively high value leading to a possible overestimation of the inundated areas, while the "MCET" algorithm identifies a final threshold at a lower value, which leads to a possible underestimation of the inundated areas, especially in emergent vegetation areas. As a sequence, the use of "avg" helps to balance between the overestimation and underestimation of inundated areas. "Alt2 (avg)" and "Alt3 (avg)" are the most accurate approaches when considering, at the same time, Camargue and Doñana marshland, with both approaches exhibiting overall Kappa over 0.81 and 0.86 for the two study areas, respectively. These areas are dominated by emergent vegetation and temporary flooded areas, and, thus, "Alt2 (avg)" and "Alt3 (avg)" could be utilized for other wetlands with similar characteristics. However, the overall Kappa is below 0.74 for both approaches when considering the Doñana complete area. The complete Doñana area is comprised of larger parts that are permanently dry. In this case, the utilization of MCET achieves a higher accuracy. Overall, "Alt2 (MCET)" and "Alt3 (MCET)" showcase a consistent overall Kappa exceeding 0.79 for all three study areas, since they avert the erroneous detection of water in dry areas. Thus, these two alternatives should be preferred for wetlands comprising large parts of dry areas.
The fact that different approaches seem to operate best in one or the other study area may be related to the models used for generating the reference maps, as well as the original field data used to build the local water presence models. In Doñana as in Camargue, the reference maps are from formulas developed to optimize water detection in the wetland types that represent very effective local conditions. Doñana marshland mostly includes seasonal marshes with relative sparse and low emergent vegetation due to the short flooding period. Tall emergent plants are rare and they are generally grazed by cattle when present [48]. Hence, the formula originally developed in Doñana does not have to perform well under dense vegetation cover. On the other hand, tall emergent plants cover large areas of semi-permanent marshes in Camargue, and the formula originally developed was specifically meant to detect water under dense and tall vegetation cover [22]. The sensors and spectral bands used to develop the original formula at each site could also influence the performance of each alternative approach tested. In Doñana, the estimation of the reference maps relied on band 5 of Landsat TM and ETM, which is similar to the SWIR-1 band of S2. The best performing alternative at this site is associated with the SWIR-1 band (e.g., Alt1 data input). In Camargue, the original formula used a combination of SWIR-2 and NIR S2 bands [22], and the three best performing alternatives use NIR in combination with SWIR-1 (e.g., Alt2 data input) or SWIR-2 (e.g., Alt3 data input) bands.
In Camargue, the accuracy of all alternative approaches decreases during winter (from December to February). This systematic error is largely due to the misclassification of the water presence under dense cover of scrubby vegetation (Salicornia marshes). This habitat, which is flooded by rainfall during this specific time of year [22], contributes to 25% of water pixels that were misclassified as dry in Figure 7. Another example is given in Figure 11, where the comparison between (a) and (b) shows that "Alt3 (OTSU)" mainly detects open water areas, and neglects the water presence under Salicornia salt marshes classified as flooded in the reference map.
During the study period, S2 satellites were passing over Camargue during late morning (between 10.10 and 10.40 CET). As a consequence, shadows, which are mainly observed in the town of Arles in the north of Camargue area, appear in S2 data. Longer shadows, which are observed from the end of autumn to early spring and mainly during winter, are misclassified as inundated areas to a larger degree with the reference map compared to the automatic thresholding alternatives. Therefore, shadow misclassification related to time of the satellite passage factor can also contribute to the reduction of agreement between the reference map and automatic thresholding approaches in the winter. Another example of disagreement is related to deep waters in a large lagoon that were misclassified as non-inundated areas by the reference map in November (14 November 2017) but not by "Alt2 (avg)" (Figure 12). This misclassification (see within the red square of Figure 12b) was attributed to the presence of waves caused by strong winds during the S2 image acquisition. Presumably, other sources of "false misclassifications" could arise from the reference maps of Doñana and Camargue because they were built from models that were not 100% accurate.
December to February). This systematic error is largely due to the misclassification of the water presence under dense cover of scrubby vegetation (Salicornia marshes). This habitat, which is flooded by rainfall during this specific time of year [22], contributes to 25% of water pixels that were misclassified as dry in Figure 7. Another example is given in Figure 11, where the comparison between (a) and (b) shows that "Alt3 (OTSU)" mainly detects open water areas, and neglects the water presence under Salicornia salt marshes classified as flooded in the reference map. During the study period, S2 satellites were passing over Camargue during late morning (between 10.10 and 10.40 CET). As a consequence, shadows, which are mainly observed in the town of Arles in the north of Camargue area, appear in S2 data. Longer shadows, which are observed from the end of autumn to early spring and mainly during winter, are misclassified as inundated areas to a larger degree with the reference map compared to the automatic thresholding alternatives. Therefore, shadow misclassification related to time of the satellite passage factor can also contribute to the reduction of agreement between the reference map and automatic thresholding approaches in the winter. Another example of disagreement is related to deep waters in a large lagoon that were misclassified as non-inundated areas by the reference map in November (14 November 2017) but not by "Alt2 (avg)" (Figure 12). This misclassification (see within the red square of Figure 12b) was attributed to the presence of waves caused by strong winds during the S2 image acquisition. Presumably, other sources of "false misclassifications" could arise from the reference maps of Doñana and Camargue because they were built from models that were not 100% accurate. Another finding is related with the rice fields within the study areas. The criteria for mapping the water-vegetation subclass (see paragraph 2.4.3) as defined for Doñana are not satisfied for Camargue for any of the dates or any of the alternative approaches. On the contrary, for Doñana, the criteria are satisfied for the dates falling within summer (e.g., 11 July, 2017 and 20 August, 2017) and the water-vegetation subclass is mainly detected in the upper east area where rice fields are located. This result complies with the results in Reference [23] and is in agreement with the growing cycle of rice, since the rice grows during the summer period in Doñana and, at the same time, the paddies are flooded [59]. The difference between Camargue and Doñana is that the thresholds examined in the criteria for detecting water-vegetation subclass can be detected in the Inp for Alt1 and Alt2 and MNDVI histograms of Doñana, but not of Camargue (see Figure 4 and Figure 5). This likely relates to the land cover synthesis of each area. In particular, for the rice paddies of Camargue, it is found that there are other land cover types (e.g., reed marsh) with similar spectral behavior (i.e., similar Inp and MNDVI values) to the rice paddies. This is evident in a much smaller degree for Doñana. Moreover, the average MNDVI value of the Camargue rice paddies is lower than the average MNDVI Another finding is related with the rice fields within the study areas. The criteria for mapping the water-vegetation subclass (see paragraph 2.4.3) as defined for Doñana are not satisfied for Camargue for any of the dates or any of the alternative approaches. On the contrary, for Doñana, the criteria are satisfied for the dates falling within summer (e.g., 11 July, 2017 and 20 August, 2017) and the water-vegetation subclass is mainly detected in the upper east area where rice fields are located. This result complies with the results in Reference [23] and is in agreement with the growing cycle of rice, since the rice grows during the summer period in Doñana and, at the same time, the paddies are flooded [59]. The difference between Camargue and Doñana is that the thresholds examined in the criteria for detecting water-vegetation subclass can be detected in the Inp for Alt1 and Alt2 and MNDVI histograms of Doñana, but not of Camargue (see Figures 4 and 5). This likely relates to the land cover synthesis of each area. In particular, for the rice paddies of Camargue, it is found that there are other land cover types (e.g., reed marsh) with similar spectral behavior (i.e., similar Inp and MNDVI values) to the rice paddies. This is evident in a much smaller degree for Doñana. Moreover, the average MNDVI value of the Camargue rice paddies is lower than the average MNDVI value of the Doñana rice paddies and, thus, the discrimination of Camargue rice paddies from other vegetated areas in the MNDVI histogram is impeded. Due to the previously mentioned reasons, emergent rice cannot be detected in Camargue from the middle of July to early September. Furthermore, by comparing Alt1, Alt2, and Alt3, it is evident that Alt2 and Alt1 histograms allow the detection of T upper (see Figures 4b and 13a), while, in the histogram of Alt3, T upper cannot be detected (see Figure 13b). As a consequence, when comparing the inundation maps presented in Figure 14, it is evident that the water-vegetation subclass, corresponding mainly to rice fields located in the upper right part of the inundation map, cannot be detected when using the "Alt3 (MCET)" approach. value of the Doñana rice paddies and, thus, the discrimination of Camargue rice paddies from other vegetated areas in the MNDVI histogram is impeded. Due to the previously mentioned reasons, emergent rice cannot be detected in Camargue from the middle of July to early September. Furthermore, by comparing Alt1, Alt2, and Alt3, it is evident that Alt2 and Alt1 histograms allow the detection of Tupper (see Figure 4b and Figure 13a), while, in the histogram of Alt3, Tupper cannot be detected (see Figure 13b). As a consequence, when comparing the inundation maps presented in Figure 14, it is evident that the water-vegetation subclass, corresponding mainly to rice fields located in the upper right part of the inundation map, cannot be detected when using the "Alt3 (MCET)" approach.  This study proves that automatic thresholding can be applied to more than one study areas and achieve high inundation mapping accuracy, without the need for simultaneous ground truth data or user's intervention. Machine learning approaches that are developed based on ground truth data derived from a specific site [13,30,33,35,37] may perform very accurately as well. However, a credible value of the Doñana rice paddies and, thus, the discrimination of Camargue rice paddies from other vegetated areas in the MNDVI histogram is impeded. Due to the previously mentioned reasons, emergent rice cannot be detected in Camargue from the middle of July to early September. Furthermore, by comparing Alt1, Alt2, and Alt3, it is evident that Alt2 and Alt1 histograms allow the detection of Tupper (see Figure 4b and Figure 13a), while, in the histogram of Alt3, Tupper cannot be detected (see Figure 13b). As a consequence, when comparing the inundation maps presented in Figure 14, it is evident that the water-vegetation subclass, corresponding mainly to rice fields located in the upper right part of the inundation map, cannot be detected when using the "Alt3 (MCET)" approach.  This study proves that automatic thresholding can be applied to more than one study areas and achieve high inundation mapping accuracy, without the need for simultaneous ground truth data or user's intervention. Machine learning approaches that are developed based on ground truth data derived from a specific site [13,30,33,35,37] may perform very accurately as well. However, a credible This study proves that automatic thresholding can be applied to more than one study areas and achieve high inundation mapping accuracy, without the need for simultaneous ground truth data or user's intervention. Machine learning approaches that are developed based on ground truth data derived from a specific site [13,30,33,35,37] may perform very accurately as well. However, a credible performance for other sites cannot be safeguarded. Several studies have identified that surface reflectance accuracy of S2 L2A products may vary for different dates [60,61]. Hence, the performance of machine learning approaches could be negatively affected for dates that there is a notable variance in the surface reflectance accuracy compared to the dates, for which the training samples were derived in order to train the classification models.

Conclusions
This study examines automatic thresholding alternative approaches for separating inundated class pixels from non-inundated class pixels by utilizing atmospherically corrected S2 data. The experimental results show that alternative approaches are able to achieve high classification accuracy for Camargue and Doñana study areas. Out of the nine alternative approaches tested, three, seven, and four approaches were "Alm Perf" for Camargue, Doñana marshland, and Doñana complete area, respectively. "Alt2 (avg)" and "Alt3 (avg)" provided "Alm Perf" results for both Camargue and Doñana marshland, while "Alt2 (MCET)" and Alt3 (MCET)" provided the most consistent results for all areas, including the Doñana complete area. Thus, "Alt2 (avg)" and "Alt3 (avg)" are suggested for wetlands extensively covered by temporary flooded and emergent vegetation areas, such as the Camargue and the Doñana marshland, while "Alt2 (MCET)" and "Alt3 (MCET)" are expected to give more consistent results for wetlands including a large portion of dry areas, such as the Doñana complete area.
Future steps could consider the exploitation of ancillary information, such as digital elevation models to improve water detection under emergent vegetation, by inferring the water presence based on detected adjacent water covered areas having similar elevation, and land cover information to correct areas erroneously classified as water covered, where water presence is not expected. Furthermore, S2 inundation maps of a site generated via the automatic thresholding alternative achieving top accuracy among other alternatives can be fused with S1 data in order to allow for inundation mapping during extended cloudy periods, based on the example of Reference [12].
Author Contributions: G.A.K. and I.M. conceived and implemented the alternative approaches in this work, which were fine-tuned based on suggestions from G.L. and B.P. G.A.K. performed an accuracy assessment of the generated Sentinel-2 inundation maps. Reference maps for Camargue were provided by G.L. and B.P. Results analysis was carried out by all authors. G.A.K. and I.M. mainly wrote the paper, with significant contribution of G.L. and B.P. to the writing of the discussion and the description of the Camargue study area. All authors reviewed, suggested improvements, and approved the manuscript.