Estimation of Water Coverage in Permanent and Temporary Shallow Lakes and Wetlands by Combining Remote Sensing Techniques and Genetic Programming: Application to the Mediterranean Basin of the Iberian Peninsula

: This work aims to validate the wide use of an algorithm developed using genetic programing (GP) techniques allowing to discern between water and non-water pixels using the near infrared band and different thresholds. A total of 34 wetlands and shallow lakes of 18 ecological types were used for validation. These include marshes, salt ponds, and saline and freshwater, temporary and permanent shallow lakes. Furthermore, based on the spectral matching between Landsat and Sentinel-2 sensors, this methodology was applied to Sentinel-2 imagery, improving the spatial and temporal resolution. When compared to other techniques, GP showed better accuracy (over 85% in most cases) and acceptable kappa values in the estimation of water pixels ( κ ≥ 0.7) in 10 of the 18 assayed ecological types evaluated with Landsat-7 and Sentinel-2 imagery. The improvements were especially achieved for temporary lakes and wetlands, where existing algorithms were scarcely reliable. This shows that GP algorithms applied to remote sensing satellite imagery can be a valuable tool to monitor water coverage in wetlands and shallow lakes where multiple factors cause a low resolution by commonly used water indices. This allows the reconstruction of hydrological series showing their hydrological behaviors during the last three decades, being useful to predict how their hydrological pattern may behave under future global change scenarios.


Introduction
The Mediterranean basin holds many wetlands of high ecological, social, and economic value. Yet, many of these important natural assets have been notably degraded or even destroyed, mainly during the 20th century [1]. The causes of the degradation of these ecosystems are multiple: uncontrolled water extractions, artificial damming of feeding watercourses, artificial water inputs, drainage structures, and aquifer overexploitation, among others [2].
The importance of Mediterranean wetlands is acknowledged in multiple multilateral environmental agreements, especially the Ramsar Convention on Wetlands [3,4], one of the most influential agreements for the conservation of wetlands globally [5]. Furthermore, several policies have been established to protect, conserve, and manage these ecosystems. The Ramsar Convention [4], as well as the Habitats Directive (92/43/EEC) and the Water Framework European Directive (2000/60/EC), and other worldwide, European, or national due to the higher sensitivity of these wavelengths to the presence of water [13,19,20], thus causing mistakes in the water-cover pixel identification.
More recently, Doña et al. [20] introduced a remote sensing methodology to estimate the changes in water coverage of lenitic (standing waters) ecosystems using Landsat-7 ETM+ images. This approach was based on an algorithm developed using genetic programing (GP) techniques that allows to discern between water and non-water pixels using the near infrared band of Landsat-7 and different thresholds. The selection of the GP algorithm was based on a higher stability of the thresholds versus other methods, where in GP, the algorithms differentiated for both playa lakes and non-playa lakes in permanent and seasonal Mediterranean saline lakes. This algorithm was first developed to be applied to shallow lakes, both permanent and temporary, of the Biosphere Reserve of La Mancha Húmeda, Spain [20]. This approach might be a useful tool to monitor the hydrological patterns of highly fluctuating seasonal and permanent lakes in semiarid areas, such as the Mediterranean basin. Consequently, the aims of this work are: (1) to assess the performance of the model and the stability of the thresholds proposed by Doña et al. [20] to estimate water coverage in lenitic ecosystems focusing on several lake and wetland types, including temporary lakes, located in semiarid regions of the Mediterranean basin of the Iberian Peninsula; and (2) to study the adaptation of this methodology to Sentinel-2 imagery improving the spatial resolution from 30 m to 10 m and the revisit time from 16 to 5 days.

Study Area
This study focuses on a set of wetlands located in the Mediterranean biogeographical region of the Iberian Peninsula ( Figure 1). According to the Ramsar Convention [4], wetlands, in a broad sense, include lakes and rivers, underground aquifers, swamps and marshes, wet grasslands, peatlands, oases, estuaries, deltas and tidal flats, mangroves and other coastal areas, coral reefs, and constructed sites, such as salt pans. A total of 34 water bodies, of 18 different wetland types [36], were studied. Many of them are included in Natura 2000 Network sites (European Habitats Directive (92/43/EEC)) and are Wetlands of International Importance for the Ramsar Convention [4]. The studied wetland types include freshwater and salt marshes, salt ponds, as well as saline or freshwater, temporary or permanent lakes, among others (Table 1). These wetlands are distributed within the Spanish regions of Andalucia, Aragón, Castilla y León, Castilla-La Mancha, Catalunya, Extremadura, and Comunitat Valenciana.  (34) This work was developed in the framework of the CLIMAWET and CLIMAWET-CONS projects. These projects aim to recognize the patterns of ecological functioning of Iberian Mediterranean wetlands and the effect of its conservation status in order to study their adaptation to climate change. This is addressed through the study, per ecological type, of its biogeochemistry under natural conditions, then related to its conservation status and management modes [37]. For all these purposes, the results of the flooded area and the hydrological pattern obtained from our method provide essential information for the conservation or recovery of these valuable systems.
Mediterranean wetlands show a wide diversity regarding its origin, geographic location, hydroperiod, water chemistry, soil or sediment characteristics, and vegetation [38]. Despite this heterogeneity, Mediterranean wetlands share a number of common features, for example their seasonality determined by the Mediterranean climate or, for many cases, the flat topography of the area. In our case, some of the studied wetlands are temporary lakes and their hydroperiods fluctuate as a result of the mixed Mediterranean-Continental semiarid climate patterns [39]. Some of them have a marginal belt of ruderal or halophytic plants, and some also have helophytic vegetation mainly in marginal areas. Furthermore, the studied water bodies are distributed along a salinity gradient, ranging from freshwater to hypersaline. Table 2 shows the main features of the studied water bodies.

Reference Data
High spatial resolution Google Earth Pro images were used as a reference dataset to validate the water area estimated by the applied methods. Boundaries of the water bodies and maximum flooded area were digitized manually using Google Earth Pro Tool. A vector file was generated for each date, and then rasterized and resampled to the spatial resolution of interest using ArcGIS Desktop, i.e., 30 m or 10 m in the case of Landsat-7 ETM+ or Sentinel-2 imagery, respectively. In situ testing of the delineated boundaries were performed within field campaigns carried out during 2017-2018. In these seasonal campaigns, a sort of ecological variables was determined, along with the study of the carbon cycle [37], and ancillary in situ determinations of wetland surface (water and vegetation cover) were also performed.

Remote Sensing Data Collection and Pre-Processing
The availability of scenes for a specific wetland site depends on several factors such as the geographical location or the weather conditions (e.g., clouds affecting images). On the other hand, climatic and meteorological parameters become essential e.g., whether a year is dryer or wetter determines the appropriate number of scenes required. For the time course of the whole year hydroperiod, the combined use of satellites with different and complementary visiting periods can guarantee a proper temporal coverage fitting with the observation of the water cover dynamics. Just in cases of massive rainfall in small shallow temporary lakes, some problems could arise because of the lack of temporal encompassing of the satellite images and the eventually rapid filling, even though this could cause a delay of just a few days in the detection timing. Both Landsat-7 ETM+ and Sentinel-2A/2B imagery were used to conduct this study. Based on the spectral matching between Landsat-7 ETM+ and Sentinel-2 sensors, the methodology to estimate water and non-water pixels was applied to both sensors. The medium-high spatial resolution of these sensors enable to study small-to-medium-sized water bodies. Moreover, Landsat images are available free of charge at the United States Geological Survey (USGS) website (http://glovis.usgs.gov), and the extensive historical archive available of these images allows studies of wetlands hydrological dynamics since 1984. Sentinel images are available in the Copernicus Open Access Hub (https://scihub.copernicus.eu/dhus/#/home), which provides complete, free and open access to Sentinels user products. The images of Sentinel-2A/2B are available since 2015 and 2017, respectively. The main features of these sensors are listed in Table 3. In brackets, lake identification number (as in Figure 1). * 2 Long, Moderate, and Short refer to temporary systems with average flooding periods of more than 9 months, between 6-8 months, and less than 6 months on average, respectively. Ephemeral hydroperiod refers to those systems with no regular flooding, i.e., not directly linked with the immediate climate but to long-term aquifer dynamics. Semi-permanent remain flooded the whole year or get dry in summer other years. * 3 L corresponds to Landsat-7 ETM+ and S to Sentinel-2A/2B imagery. * n.a. non-available. † This is a naturally temporary lake that has been morphologically altered. A total set of 53 images, of dates with reference data available, were selected for this work (27 Landsat-7 ETM+ and 26 Sentinel-2A/2B). Satellite imagery dates were selected concurrent to reference data availability (Tables S1 and S2). Only satellite images close in time with the reference data were considered to reduce the errors associated to water cover changes. Hence, coincident images (same day) were considered for lakes/wetlands with short hydroperiod and a temporal window of a week was considered for lakes/wetlands with permanent hydroperiod. In the case of Landsat-7 ETM+, images corresponding to Level 2 product were downloaded. This product provides information on surface reflectance, already corrected atmospherically by means of the radiative transfer code 6S (sr) following Vermonte et al. [40]. The Landsat paths and rows covering all locations in our study area are: 198/32-33, 199/35, 200/31-34, 201/32-33, and 202/31,34. The scan line corrector (SLC) failure affecting the Landsat-7 ETM+ images was corrected following the methodology proposed by Scaramuzza et al. [41]. In the case of the Sentinel-2 imagery, the atmospherically corrected, Level-2A surface reflectance product, was also used. Sentinel bands 11 and 12 were resampled to 10 m pixel size, using the technique of the next neighbor, which is the least changing pixel value original. After this pre-processing, all images were ready for the application of the method designed to estimate the water area.

Water Mapping Method
As mentioned above, in this work we explored the robustness of the methodology developed by Doña et al. [20], to detect and monitor water cover within a sort of Mediterranean shallow lakes and wetlands. However, the novelty of this paper goes further, into the intercomparison of the performance with a wide range of other indicators and water indices used in the literature, as well as for the variety of wetlands considered throughout the Mediterranean region for this study. This approach is based on an algorithm developed using genetic programing (GP) techniques that allows to discern between water and non-water pixels using the near infrared band of Landsat-7 ETM+ and different thresholds, both for playa and non-playa lakes (Equation (1)). Based on its accuracy assessment by intercomparison of the performance with a wide range of other indicators and indices used in the literature, such as water and vegetation indexes, single band thresholds, supervised and unsupervised classification methods, and machine learning methods such as artificial neural networks and support vector machine, the GP algorithm was selected as the most appropriate method.
The estimation of water pixels across a set of 13 shallow lakes, most of them temporary, showed a kappa value of 0.98 and a low error of 2%, both in training and validation. The selection of the GP algorithm was precisely based on a higher stability of the thresholds for the assayed approach versus other methods. Doña et al. [20] reported the need to modify the threshold for the GP algorithm only under the presence of playas in the studied water bodies. Values < −0.05 for water pixels in playa lakes, and <0.1 in the case of lakes without playas, were obtained in that work.
In the present study, the performance of the GP algorithm was also compared to the results obtained through water indices techniques in a comparative study. The evaluated indices were the Normalized Difference Water Index (NDWI) [32]; the Automated Water Extraction Index (AWEI) proposed by Feyisa et al. [35], in its two versions, AWEInsh for areas not affected by shadows, as well as AWEIsh for areas with presence of shadows or built-up surfaces; the Water Ratio Index (WRI) [42]; and the Water Turbidity Index (WTI) [43]. A confusion matrix was obtained for each of the applied methods. This matrix is a standard tool for evaluating statistical models by means of the kappa coefficient (κ). This index varies from −1 to 1, where values higher than 0.4 are considered acceptable [44]. The confusion matrix is also used to generate the user commission and omission error, as well as the accuracy. The commission error is the percentage of pixels wrongly assigned to a certain class by the classifier, while user accuracy is the probability that a pixel assigned to a class by the classifier correctly corresponds to that class. The omission error is the percentage of pixels that actually belong to the class but were incorrectly classified, and the accuracy of the producer is the probability that the classifier has correctly assigned a class given by the reference data.

Landsat-7 ETM+ Imagery
Due to the spatial resolution of Landsat-7 ETM+, the presence of clouds in the scenes, and the availability of coincident images with the reference data, only 11 of the 18 studied wetland types were analyzed using Landsat-7 ETM+ imagery (Table 1). Table S1 shows the acquisition dates for the Landsat-7 ETM+ images used in this work. The total number of pixels classified as covered by water by each of the different methods studied showed a wide range of results in comparison to all the reference data set, composed of a total of 13,583 water pixels ( Figure 2). In order to determine the optimal threshold for each water index tested, multiple thresholds to differentiate between water and non-water pixels were considered. The NDWI and WRI indices showed the best results with the thresholds proposed by the original authors to assign water pixel values, i.e., values > 0 in the case of the NDWI and values > 1 for the WRI [32,42]. The threshold obtained for the AWEIsh has a good agreement with results obtained by the authors [35], showing values > 0 to assign water pixels. In the case of the AWEInsh, the best results were obtained with a threshold of −0.45 and 0.2 for the WTI. Overall, the best performance was found when using the GP0.1 method proposed by Doña et al. [20] and the AWEInsh method proposed by Feyisa et al. [35], showing kappa values of 0.6 in both cases, and overall accuracy values of 85% and 90%, respectively ( Figure 2). As shown in Figure 2, both methods presented a slight overestimation in the number of water pixels as compared with the reference data. The GP-0.05 and WTI also showed an overestimation in the number of water pixels. However, in this case, both the kappa coefficient and overall accuracy were much lower, with values of 0.4 and 75%, and 0.2 and 70%, for GP-0.05 and WTI, respectively. Figure 2 also shows that the NDWI, AWEIsh, and WRI results showed the opposite trend than the other methods, i.e., these methods showed an underestimation in the number of water pixels as compared with the reference data. In this case, the AWEIsh showed kappa and overall accuracy values of 0.5 and 85% followed by the WRI and NDWI with kappa and overall accuracy values of 0.4 and 80%, and 0.3 and 75%, respectively.  Figure 3 shows the results, sorted by lake/wetland type, of both the kappa coefficient and the overall relative accuracy (OA), including all the methods tested. As shown, GP-0.05, GP0.1, and AWEInsh showed the best performance in the estimation of water pixels for the studied lake/wetland types (with mostly k > 0.7 and OA > 90%) e.g., (i) GP's for type 1.3.2.5.1C/R-Temporary shallow hypo-mesosaline lakes; (ii) GP0.1 and AWEInsh for type 1.3.2.6.2-Non-saline shallow ponds and wetlands with alkaline waters, temporary and; (iii) AWEInsh for type 1.3.2.7.2-Non-saline shallow ponds and wetlands (morphostructural origin) with low alkalinity waters,. Furthermore, it is noteworthy that AWEInsh shows, in general, acceptable results for non-saline lake/wetland types, with k > 0. 6  With regards to the GP methods, the GP0.1 behaved better than the GP-0.05, showing kappa coefficient values above 0.8 and OA values close to 100% in 4 of the 12 assayed lake/wetland types. For GP-0.05, kappa showed values above 0.8 only for 3 types. Aside from non-saline lake/wetland type, the GP in its two versions showed acceptable results for temporal and permanent saline lake/wetland types (k > 0. 6  In general, for saline playa lakes, GP-0.05 is expected to give better results than GP0.1. Despite these differences as a function of the threshold, both GP methods yielded better results than the rest of the methods for hypo-mesosaline lakes of type 1.3.2.5.1. (both well-preserved, C, and restored, R).  Table 1, for the Landsat-7 imagery. In blue is indicated the number of reference water pixels and the best results for each type are marked in red.
Similar results (Figure 3) were obtained for the Mediterranean coastal lagoons (type 2.1.3.4.2) and the naturalized or restored salt pans (type 2.1.3.5.1). Although all methods, except GP-0.05 and WTI, presented acceptable kappa values in these cases (k > 0.4), the GP0.1 and the AWEIsh showed higher kappa and OA values (k > 0.6 and OA > 70%) for both types. However, low kappa values (k < 0.4), indicative of poor performance, were found for GP0.1 and AWEIsh for type 1.3.2.5.3, corresponding to temporary soda lakes, whereas the AWEIsh additionally showed kappa values lower than 0.4 for the type 1.3.2.1.1 corresponding to fluvial ponds and wetlands in middle-lower reaches on floodplains. Although GP0.1 and AWEInsh presented kappa values close to 0.6 for type 1.3.2.1.1, the other methods showed kappa values close to 0. These results could be due to the effects of the strips created by the SLC failure in the original ETM+ images, which affected the images used for the studied lake. Regarding the type 1.3.2.5.3 (soda lakes), among the studied lakes this includes exclusively Laguna de Caballo Alba, a very shallow lake with typical depth values of just a few centimeters. In this case, the poor performance of all the tested methods is likely linked to the limitation in discerning between water and non-water pixels for extremely shallow waters (Figure 4a). The classification outputs of AWEInsh, AWEIsh and WRI (Figure 4a.4-a.6) showed that these methods were capable of estimating the presence of water in Laguna Caballo de Alba, albeit with a considerable underestimation in the classification of water pixels in comparison with the reference image ( Figure 4a). Meanwhile, GP-0.05, GP0.1, NDWI, and WTI could not discern between water and non-water pixels (Figure 4a.1-a.3 and a.7) in this lake type (1.3.2.5.3). Typically accepted errors in wetland water coverage assessment studies at a single site range 10-20% [5]. According to this assessment, Table 4 lists the methods showing acceptable mean uncertainties in the estimation of water pixels for each lake/wetland type using Landsat 7 imagery. Total error refers to the sum of the commission and omission errors. Following this criterion, eight lake/wetland types presented total errors lower than 20%. The GP method showed acceptable total error values ranging from 3% to 17% for seven of these eight types. Results obtained for each GP version showed acceptable total error values for four and six types, for GP-0.05 and GP0.1, respectively (both GP-0.05 and GP0.  It is noteworthy that the GP method, including both versions, showed total error values above 20% only in one of the eight types. This occurred for the type 1.3.2.7.2 (temporary non-saline shallow ponds and wetlands with low alkalinity waters, morphostructural origin), which includes, among others, Laguna Grande de la Puebla de Beleña. This poor performance was caused by an overestimation of water pixels for this water body, mainly when the lake was dry. This might be due to the large amount of organic matter in its sediments. This characteristic causes misinterpretation by the algorithm as consequence of the dark color of the lake bottom when it is dry. In this case, AWEInsh was the sole method capable of detecting the dry condition of the lake. However, AWEInsh showed a total error above 20% for types 1.3.2.1.3 (permanent lakes formed by river damming in upper course), represented by Laguna del Arquillo (Figure 4), and 1.3.2.5.1.C (temporary hypo-mesosaline shallow lakes, well preserved). In this case, AWEInsh underestimated the water pixels for Laguna del Arquillo (Figure 4b), whereas this method was unable to detect total dryness in Laguna Salada de Zorrilla.

Sentinel-2 Imagery
As a result of the availability of coincident reference data and the presence of clouds, a total of 10 lake/wetland types could be finally tested using Sentinel-2 imagery (Table 1). In this case, 3 of the 10 lake/wetland types coincide with the types evaluated with Landsat-7 ETM+ imagery (1.3.2.1.3, 1.3.2.7.2, and 2.3.4.2 corresponding to fluvial permanent ponds and wetlands originated by natural damming, on upper reaches, non-saline shallow ponds and wetlands-morphostructural origin-with low alkalinity waters, temporary, and Mediterranean coastal lagoons, respectively). Table S2 shows the acquisition dates for the Sentinel-2A/2B images used in this work. Results obtained for GP0.1 and WRI to classify water pixels showed similar results as those obtained with Landsat-7 ETM+ imagery ( Figure 2) for these methods, with kappa coefficient and overall accuracy values of 0.6-85% and 0.4-85%, respectively. In the case of all other tested methods, kappa coefficient and overall accuracy values showed, in general, lower values than those obtained with Landsat-7 ETM+ imagery. In this case, results showed a clear overestimation of the number of water pixels estimated by GP-0.05, AWEInsh and WTI methods, in comparison with the reference data composed of a total of 76,511 water pixels. As in the Landsat-7 ETM+ results, the GP0.1 developed by Doña et al. [20] yielded the best kappa and overall accuracy values. Furthermore, these results agree with those obtained with Landsat-7 ETM+ imagery, as mentioned above, where the same kappa value of 0.6 and overall accuracy values of 85% were achieved for GP0.1. Contrastingly, for the GP-0.05, AWEInsh, AWEIsh and WRI methods, kappa and overall accuracy showed lower values than those obtained by applying these methodologies to Landsat-7 ETM+ imagery. Kappa and overall accuracy values decreased from 0.4 to 0.3 and from 75% to 60%, respectively, for GP-0.05, from 0.5 to 0.4 and from 85% to 80% in the case of AWEIsh, while they dropped from 0.6 to 0.5 and from 90% to 80% for AWEInsh. In contrast, the NDWI increased its kappa and overall accuracy values from 0.3 to 0.5 and from 75% to 90%, respectively, when applying this method to Sentinel-2 imagery compared to Landsat-7 ETM+. Figure 5 shows the results, in terms of both the kappa coefficient and the overall accuracy for the different lake/wetland types. In general, GP0.1 showed the best results compared to the other tested methods. GP0.1 showed kappa coefficient and overall accuracy values above 0.6 and 80% for 6 of the 10 evaluated types. As shown in Figure 5, GP0.1 and AWEInsh showed the best performance in the estimation of water pixels for shallow lake/wetland types (k > 0.6 and OA > 80%), e.g., GP0.1 and AWEInsh for type 1.3.2.6.1-non-saline shallow ponds and wetlands with alkaline waters, permanent, and GP0.1 for type 1.3.2.7.2-non-saline shallow ponds and wetlands (morphostructural origin) with low alkalinity waters, temporary. These results are in agreement with those obtained with Landsat-7 ETM+ imagery. Moreover, the GP0.1 showed acceptable results for temporal/permanent and saline/non-saline lake/wetland types (k > 0. 6 Table 1, for the Sentinel-2 imagery. In blue is indicated the number of reference water pixels and the best results for each lake/wetland type are marked in red. The AWEIsh, NDWI, and WRI water indices showed similar results among them with acceptable kappa coefficient values (k > 0.5) and overall accuracy values over 70% for types 1.  (Figure 5). In these cases, both GP methods, as well as AWEInsh and WTI, showed an overestimation of water pixels for all studied water bodies. The opposite trend was found for these types for the NDWI, AWEIsh, and WRI methods. On the other hand, all the methods, except NDWI, presented an overestimation of water pixels in Mallada de la Mata del Fang (l'Albufera) (intradunal ponds, type 2.1.3.1.2.6). In this case, the NDWI was unable to discern between water and non-water pixels, as it got that the water body was permanently dry. Figure 6b shows the thematic maps obtained for Filtre Biologic l'Embut (type 2.1.3.4.1.) where the differences among the results of the application of the different methods tested with Sentinel 2 imagery are displayed. Although the thematic maps, at least for the case of NDWI, AWEInsh, AWEIsh, and WRI (Figure 6b.3-b.6), are apparently good, the total error obtained in the classification of water pixels using these methods was around 30%. Furthermore, the type 1.3.2.7.2 (temporary non-saline shallow lakes and wetlands with low alkalinity) showed unacceptable kappa values for all the assayed methodologies, using Sentinel imagery, except for the GP0.1. This comes mainly from an overestimation of water pixels in Charca del Monte and Laguna de Campillo. For the latter, the NDWI, AWEIsh, and WRI could not even discern between water and non-water pixels, and they described these lakes as dry (Figure 6a.3, a.5 and a.6) when they actually presented some flooding (Figure 6a).  Table 5 shows the methods showing acceptable mean uncertainties (lower than 20%) in the estimation of water pixels for each lake/wetland type when applied to Sentinel-2A/2B imagery. Total error refers to the sum of the commission and omission errors. Only 4 types, of up to 10 assayed, presented total errors lower or equal than 20% ( Table 5). The GP0.1 method showed acceptable total error values ranging from 11% to 20% for these four types, followed by the AWEInsh and the NDWI for which the total error values range from 13% to 20% for two types. In the case of AWEIsh and WRI, the results showed acceptable total error values for two and one types, respectively. Contrasting to the Landsat-7 ETM+ results, the Mediterranean coastal lagoons (2.1.3.4.2), that includes L'Encanyissada, Senillar de Moraira, and Albufera Nueva de Adra, showed unacceptable values in the classification of water pixels, with total error values above 20% in the classification results when all the different methods were applied to Sentinel 2 imagery. This result could be due to the low kappa values shown by Senillar de Moraira, which were lower than 0.4 for all the assayed methods, mainly due to the small flooded surface corresponding only to 10 pixels, most of them mixed (water, vegetation, and soil). The analysis of the results for 1.3.2.1.3 (permanent lakes formed by river damming in upper course) is of particular interest due to the availability of coincident Landsat-7 ETM+ and Sentinel-2 A/B images for this type. It is noteworthy that the application of the different methodologies for the type 1.3.2.1.3 responded in a similar way when using both Landsat-7 ETM+ and Sentinel-2A/2B imagery (Figure 7). GP0.1, GP-0.05, and AWEInsh were the only methods capable of estimating the water pixels for this type 1.3.2.1.3 with acceptable errors (<20%). Very similar results were obtained for GP and AWEInsh, with kappa values of 0.9 in the estimation of water pixels. Water cover was underestimated in this case by both GP0.1 and AWEInsh approaches, whereas the opposite trend, overestimation, was observed for GP-0.05 (Figure 7). This was not the case for the type 1.3.2.7.2., though for this type, the studied lakes as determined by the availability of images and reference data were different i.e., those studied with Landsat-7 ETM+ imagery were Laguna de Talayuelas and Laguna Grande de la Puebla de Beleña, while those studied with Sentinel-2 imagery were Laguna de Campillo, Laguna de Palancoso and Charca del Monte.

Discussion
In this study, we widely tested the robustness of the methodology proposed by Doña et al. [20] a specific algorithm developed by genetic programming (GP) which allows to discern between water and non-water pixels by use of NIR band of Landsat-7 ETM+ sensor. In general, the GP was shown to assess flooding accurately, particularly in temporal shallow lakes/wetlands. These results are consistent with those published by Doña et al. [20], as GP was developed from shallow lakes, most of them temporary. Furthermore, the GP was capable of distinguishing better than other methods between water and non-water pixels for some saline/non-saline and temporal/permanent lake/wetland types. GP showed acceptable kappa coefficient values (K > 0.6) and overall accuracy values above 90% for 7 of the 11 lake/wetland types evaluated with Landsat-7 ETM+ imagery (types 1.  [20] where they observed different thresholds between non-playa and playa lakes. The comparison among the GP methodology with several standard water indices found in the literature, such as NDWI [32], AWEI [35], WRI [42], and WTI [43], revealed that GP and AWEI methods estimate the surface water in the studied lake/wetland types more accurately than the other tested methods, especially for temporary shallow lakes. GP was originally developed and tested across several water bodies including temporary shallow lakes [20], whereas AWEI was tested by its authors under a wide range of environmental conditions and water body types, such as small freshwater reservoirs, large lakes, harbors, and the sea [35]. In our study, AWEI presented better results for permanent lake/wetland types than for temporary lake/wetland types, which could possibly be explained by the fact that it was developed and tested with permanent water bodies [45][46][47][48]. In general, AWEInsh behave better than AWEIsh, mainly due to the absence of shadows in the scenes [35]. However, AWEIsh showed better accuracy results in the estimation of water pixels (k = 0.8 and OA = 92%) than those of AWEInsh in Palmones marshes site. This result agrees with those obtained by Fisher et al. [46], where AWEIsh achieved greater accuracy in the estimation of water pixels for images with urban areas, as is the case of Palmones marshes site.
NDWI, WRI, and WTI indices assessed surface water coverage with lower accuracy than GP and AWEI methods [12,20,35,46]. The water cover estimation errors by the former were basically linked to the spectral characteristics of the water bodies themselves. As these water indices were developed from permanent waters, the spectral features of this water bodies are close to the typical spectral reflectance for clear water, with low concentration of sediments compared to the high rates characteristic of the seasonal shallow water bodies [13].
Knowing that typically accepted errors in wetland water coverage assessment studies at a single site range 10-20% [5], among all methods tested, GP was the method that showed the best accuracy in the estimation of water pixels, showing good results for 7 of the 11 lake/wetland types evaluated with Landsat-7 ETM+ imagery and for 4 of the 10 lake/wetland types evaluated with Sentinel-2 imagery. For all the applied methods, results showed some discrepancies (k < 0.4) in the estimation of water cover in lakes/wetlands with presence of high coverage of helophytic vegetation and an irregular and highly fluctuating distribution of flooded areas along the water bodies. Recently, Pena-Regueiro et al. [46] found that the best performing index for four representative coastal wetlands with this characteristics was the NDWI with a threshold value of −0.3, which showed overall accuracy and kappa coefficient results above 75% and 0.5 for all the wetlands studied, respectively.
The spatial resolution of Landsat-7 ETM+, 30 m, facilitates the mixture of land-cover types, this effect being more significant in edge pixels, resulting in difficulties for mapping flooded areas in heterogeneous wetlands [12,20,46]. Furthermore, as was mentioned in the results section, the SLC failure of Landsat-7 imageries caused errors in the estimation of the total water area of some lakes and lakes/wetlands presenting irregular shapes could be more vulnerable to this problem. However, Sentinel-2, with higher spatial resolution (10 m) than Landsat 7-ETM+, is capable of defining with more accuracy the boundary of water surfaces reducing the number of pixels with a mixture of land covers such as soil, vegetation, or water. In other cases, such as the application of the different methodologies for the type 1.3.2.1.3-permanent lakes formed by river damming in upper course, the performance using both Landsat-7 ETM+ and Sentinel-2A/2B imagery was similar, as published by Lefebvre et al. [14], who identified similar threshold values for Landsat and Sentinel-2 imagery in several land cover types in the Camargue wetlands (southern France).
GP was tested for several lake/wetland types, showing that chemical and biophysical features in the water bodies could affect the performance [20]. Although better results in the discrimination between water and non-water pixels were obtained by the GP algorithm in comparison with other water indices in this work, further evaluations and refinements of the methodology are required to widen its effectiveness under a variety of water and environmental conditions. More detailed studies are needed for the intradunal ponds, the Mediterranean microtidal swampy plains, and the Mediterranean marshes not connected to the sea, where the thresholds proposed for GP were unable to accurately distinguish between water and non-water pixels, mainly due to the presence of high coverage of helophytic vegetation and an irregular and highly fluctuating distribution of flooded areas along the water bodies.
Mapping flooded areas is essential, since wetland biodiversity reacts to flooding regimes [49][50][51], both in terms of water cover area and hydrological cycles. Wetland's biogeochemistry, including the carbon cycle and their exchanges of greenhouse gases (GHG) with the atmosphere, are also intimately linked to the water regime [52,53]. Overall, the hydroperiod is a crucial indicator for ecological trends in ecology features, functions and services of Mediterranean wetlands [5]. The use of Landsat imagery is especially attractive due to the extensive historical archive available that allows for the reconstruction of the hydrological patterns of wetlands over the last three decades. This might be useful for describing the hydrological patterns of these lakes/wetlands related to the annual meteorological pattern, among other reasons, and this, in turn, allows to forecast their behavior in future climate scenarios. Furthermore, the spectral band similarity with Sentinel-2 showed the continuity to the GP algorithms designed for Landsat imagery. Moreover, applying this methodology to Sentinel-2 imagery offers the opportunity to continue the monitoring of hydrological patterns, but improving the revisit cycle from 16 to 5 days and the spatial resolution from 30 m to 10 m. The higher spatial and temporal resolution gained with Sentinel-2 is especially important to monitor the hydrological patterns of patchy and highly fluctuating seasonal and permanent lakes and wetlands in semiarid areas, such as the Iberian Mediterranean basin. Note that a revisit frequency of five days might not be sufficient to capture the water changes in extremely dynamic water bodies, but it can fairly accomplish the monitoring requirements for the majority of lakes/wetlands types, even those temporary, in the Mediterranean Basin of the Iberian Peninsula. This might strengthen management strategy-linked lake conservation, and help to accomplish the goals of both the European Water Framework Directive (2000/60/EC) and the Habitats Directive(92/43/EEC), but also could be useful for the improvement of its role as carbon sinks under a climate change scenario [52,53] where the role of flooding, as well as that of macrophytes [49][50][51], would be of great relevance.

Conclusions
An approach for water mapping based on an algorithm developed using genetic programing (GP) developed by Doña et al. [20] has been evaluated in this study, applied to Landat-7 ETM+ and Sentinel-2A/2B imagery. Results show the potential of this method to be used for a wide variety of lakes and wetlands within the Mediterranean basin, improving traditional water mapping methods using water indices, especially for temporary shallow lakes and wetlands. The assessment of the GP approach across a set of 34 water bodies of 18 lake/wetland types in the Mediterranean biogeographical region of the Iberian Peninsula resulted in kappa values above 0.7 and total error values lower than 20% in the estimation of water pixels for 10 of the 18 types evaluated with Landsat-7 ETM+ and Sentinel-2 imagery. We demonstrated that GP methods improve very much the water cover determination by remote sensing methods of shallow temporary waterbodies in comparison to other methods that have usually been developed for permanent waterbodies. Another important outcome of this work is the adaptation of the GP methodology to Sentinel-2 imagery, reproducing results obtained with Landsat images. The reduction of the satellite revisit time to 5 days and the pixel size to 10 m multiplies the application possibilities of the exposed methodology. For this purpose, the GP methodology will need to be further refined for specific wetland types, such as Mediterranean coastal marshes and humid dune slacks. Furthermore, the GP approach could be transferable to other already existing and upcoming, satellite sensors containing the NIR spectral band. GP algorithms applied to remote sensing satellite imagery could serve as a valuable tool to monitor water coverage in wetlands and lakes of different ecological types, which is very promising for the future water resources management. This brings also a unique opportunity to improve our understanding of the temporal trends in seasonal lakes and their dependence on meteorological patterns, in order to establish sustainable water resources management plans that preserve the ecological health and the ecological services of wetlands and shallow lakes in semiarid basins, such as the Iberian Mediterranean basin.
Supplementary Materials: The following are available online at https://www.mdpi.com/2072-429 2/13/4/652/s1, Table S1: Landsat-7 ETM+ and reference data dates used to carry out the study, and number of water bodies for each date., Table S2: Sentinel-2A/2B dates used to carry out the study and number of water bodies for each date.  Data Availability Statement: Data are available at: https://www.dropbox.com/home/Remote%20 Sensing%202021_GPVal upon request to: carolina.dona@uv.es.

Conflicts of Interest:
The authors declare no conflict of interest.