*4.1. Supervised Classification*

Supervised machine learning algorithms establish relationships between input variables and target prediction [37–39]. We compared two popular machine learning algorithms, namely random forest (RF) and support vector machine (SVM), as supervised classifiers of surface water inundation, as predicted by the multisensory composite of Sentinel-1, Sentinel-2, NDVI, and NDWI (Table 1). These algorithms, which are available as functions within the GEE platform, were trained using the training subset, where the ground truth data served as a binary categorical response variable (0 = not an inundated wetland, 1 = inundated wetland). The trained algorithms were tested using the test subset, which was withheld from model fitting, and the model with the best performance was selected for the classification of the multisensory composite to identify the water bodies across the study area. The best performance was identified through accuracy assessment using statistical coefficients.

RF is an advanced version of a decision tree algorithm. Decision tree algorithms are robust predictive machine learning models that utilize a tree structure to establish relationships between inputs and outcomes. A tree structure mirrors how a tree starts at a wide trunk and splits into smaller branches as it is developed upward. Likewise, a decision tree learner uses a structure of branching decisions that lead examples into a final predicted class value. RF improves decision trees by combining bootstrap aggregation with random feature selection to add additional diversity to the model [40]. Further, though a decision tree is constructed on a whole dataset using all the features of interest, RF randomly selects observations and specific features to create multiple decision trees, and then averages the results to make predictions, which results in a more robust model [41]. The hyperparameters of the RF, including the number of trees, min leaf population, and bag fraction, were determined through a trial–error procedure in which we added the values gradually to obtain the least error values in the training data prediction outcome. The optimum hyperparameters of the RF model in this study are presented in Table 2.

**Table 2.** The optimum hyperparameters of the SVM and RF algorithms used for surface water classification.


The SVM classification tool uses machine learning theory to maximize predictive accuracy while automatically avoiding over-fitting the data [42]. SVM can be defined as systems that use the hypothesis space of linear functions in a high dimensional feature space, trained with a learning algorithm from the optimization theory [43]. SVM can be imagined as a surface that creates a boundary between plotted points in a multidimensional space representing their feature values. An SVM's goal is to create a flat border, called a hyperplane, which divides the space to develop relatively homogeneous partitions on either side. We adjusted the hyperparameters needed for SVM through a trial–error procedure to identify the optimum structure of the SVM model (Table 2).

#### *4.2. Identification of Wetland Surface Water*

We used the 2017 testing data to estimate wetland surface water via the trained algorithm without refitting the model. As we mentioned before, the JRC product was used to exclude the permanent water pixels in order to identify surface water in wetlands across the test site. The remaining water pixels were stratified based on the presence of emergent vegetation, allowing us to determine the accuracy of detecting surface water in vegetated vs. non-vegetated wetlands, where vegetation status was inferred from NDVI values. The Jenks natural breaks optimization method was used to classify the wetlands into three low, medium, and high NDVI clusters (Table 3). The Jenks method is a data clustering technique designed to determine the best combination of values into different classes. This is performed by attempting to minimize the variance within classes, and maximize the variance between classes. NDVI values below zero typically represent open water [44], and increase with increasing vegetation cover until they saturate for high vegetation closed canopies [45].


**Table 3.** NDVI cut-off values for classifying vegetation status in the identified wetlands.

#### *4.3. Accuracy Assessment*

The SVM and RF algorithms used to classify the multisensor composite were evaluated by constructing a confusion matrix for each model. The accuracy assessment was performed on the test subset in which the predictions and the ground truth data were compared using statistical coefficients (Equations (1–4)). Accuracy assessment was also carried out for the next year (2017) to evaluate the generalizability of the optimum model. The ground truth points for the year 2017 were only used for testing the model. Additionally, we performed separate accuracy assessments for small vegetated and small non-vegetated wetlands. To assess the accuracy of the methodology on small and highly vegetated wetlands, we provided a test set of 679 random points from small wetlands in the study area, with an additional 763 points from non-wetland classes from the 2016 aerial survey inventory. These were novel points that were not used as part of model training. These points came from wetlands with areas that ranged between 10 to 850 m2, and the presence of emergent vegetation ranged from 40 to 100%. We used the same procedure for small (10 to 850 m2) non-vegetated wetlands by providing a test set of 1680 points (1311 wetland and 369 non-wetland classes).

$$Accuracy = \frac{TP + TN}{TP + TN + FP + FN} \tag{1}$$

$$Specificity = \frac{TN}{TN + FP} \tag{2}$$

$$Sensitivity = \frac{TP}{TP + FN} \tag{3}$$

In the equations above, *N* indicates the total number of observations; *n* denotes the number of accurately classified wetland and non-wetland pixels; *TP*, *TN*, *FP*, and *FN* refer to true positive, true negative, false positive, and false negative, respectively.

$$Kappa = \frac{P\_o - P\_\varepsilon}{1 - P\_\varepsilon} \tag{4}$$

where *po* is the relative observed agreement, and *pe* is the hypothetical probability of chance agreement:

$$P\_o = \frac{TP + TN}{n} \tag{5}$$

and

$$P\_{\varepsilon} = \frac{1}{\sqrt{N}}((TP + FN)(TP + FP) + (FP + TN)(FN + TN))\tag{6}$$

### **5. Results**

The trained SVM and RF algorithms were used to classify multisensor composites for the years 2016 and 2017. The accuracy assessment showed SVM and RF models yielded favorable results across the testing data. However, the RF outperformed the SVM in both 2016 and 2017 testing data. Therefore, the RF model was selected as the optimum model for wetland inundation mapping. The overall testing data accuracy for the SVM and RF model for the year 2016 was 0.88 and 0.95 (Table 4), and for the year 2017 was 0.88 and 0.94, respectively (Table 5). A summary of accuracy assessment using overall accuracy, Kappa, Sensitivity, and specificity for the years 2016 and 2017 is shown in Tables 4 and 5, respectively.

**Table 4.** Accuracy assessment of supervised classification of wetland surface water (WSW) vs. other classes (OC) for the year 2016 on the test data subset that was withheld from model calibration.


**Table 5.** Accuracy assessment of supervised classification of wetland surface water (WSW) vs. other classes (OC) for 2017 on the test data subset that was withheld from model calibration.


We mapped wetland surface water across the study area for the years 2016 and 2017. Figures 6 and 7 show the identified wetlands across the study area for the years 2016 and 2017, respectively. The spatial resolution of the final maps is 10 m. The local wetland inundation in the study area can also be extracted based on the results. For example, a portion of the study area is magnified in Figure 8. A visual comparison between the aerial survey (ground truth data) and wetland surface water map (based on RF classifier) in Figure 9 shows that surface water in wetlands was mapped with acceptable accuracy (overall accuracy: 0.95; Kappa: 0.9). As we mentioned before, we also tested our algorithm in the identification of surface water in small vegetated and small non-vegetated wetlands. We compared the results with NDWI and Landsat-derived JRC surface water products (Tables 6 and 7). The results showed higher accuracy in RF as the optimum model (overall accuracy 0.76) compared to JRC (overall accuracy 0.60) and NDWI (overall accuracy 0.62) in surface water detection in small and highly vegetated wetlands (Table 6). The RF (overall accuracy 0.81) also outperformed the NDWI (overall accuracy 0.44) and JRC (overall accuracy 0.41) in small non-vegetated wetlands (Table 7).

**Figure 6.** Spatial distribution in inundated wetlands (**A**), spatial distribution of surface water (**B**) in August 2016.

**Figure 7.** Spatial distribution of emergent vegetation in inundated wetlands (**A**), spatial distribution of surface water (**B**) in August 2017.

**Figure 8.** Visual comparison of wetland inundation maps between predicted (**A**) and observed wetlands in aerial surveys (**B**) in a portion of the study area. The accuracy assessment for small vegetated wetlands and small non-vegetated wetlands are presented in Tables 6 and 7, respectively.

**Table 6.** Accuracy assessment for detection of wetland surface water (WSW) vs. other classes (OC) in wetlands that are both small (<850 m2) and highly vegetated wetlands (vegetation > 40%) from the year 2016 on the test data subset that was withheld from model calibration.


**Table 7.** Accuracy assessment for detection of wetland surface water (WSW) vs. other classes in small (<850 m2) and non-vegetated (vegetation < 40%) wetlands for 2016.


Figures 6 and 7 show the identified wetlands after excluding the permanent wet pixels for the years 2016 and 2017, respectively. The presence of emergent vegetation within the identified wetlands, as indicated by NDVI, for both years is also shown in Figures 6 and 7.

**Figure 9.** The result of wetland inundation maps for the year 2017 was obtained from the multisensor composite classification using random forest. The image shows the identified wetland. The background image is Sentinel-1 VV. The black arrows show some examples of those small wetlands that were not detected in the JRC (**A**); frequency of surface water occurrence from the year 1984 to 2019 obtained from the Landsat-derived JRC products (**B**); surface water visualization using the Sentinel-2 derived NDWI (**C**); water extent for the year 2017 derived from the Landsat-derived JRC product (**D**).

#### **6. Discussion**

This study developed an automated workflow within the GEE platform for mapping wetland surface water for 2016 and 2017 by applying the RF classifier to a combination of Sentinel-1, Sentinel-2 band data, and spectral reflectance indices derived from Sentinel-2. The results were evaluated using statistical coefficients and visual comparison with ground truth data, as well as results from Landsat-derived surface water products. The inundation of relatively large and deep water bodies can be identified in most existing remote sensing products. However, mapping wetland surface water in the PPR region is challenging due to two main reasons: (1) most PPR wetlands are very small and are highly sensitive to climate variability; and (2) the wetlands can be dry or wet, and they can contain different species of vegetation that can mask surface water. Therefore, these wetlands have complex spectral characteristics that complicate the detection of surface water extent from satellite sensors. Our approach also provides information regarding emergent vegetation within those wetlands. This is important because emergent vegetation provides shelter and food for aquatic vertebrates, such as waterfowl communities [46]. Our method can also detect water below those vegetation canopies; water that would otherwise be excluded from

habitat maps. We also provide an open-science algorithm in GEE for repeating these estimates, which can form the basis of long-term wetland surface water monitoring in the PPR.

A typical approach for mapping wetlands uses passive remote sensing that relies on water's optical properties, which differ from other land use types [47–49]. For instance, water quickly absorbs electromagnetic radiation, and more rapidly attenuates longer wavelengths than shorter ones [50–52]. However, the application of optical sensors in identifying PPR wetlands is limited, since both water depth and mixed pixels can change the water spectral signature [50–52]. Moreover, organic carbon compounds, water turbidity, chlorophyll content, and suspended materials can also add variation to water spectral properties. We addressed this issue by integrating the high-resolution bands of optical and radar sensors. Figure 9 shows a visual comparison of surface water derived from different remote sensing data. The figure shows that many small wetlands were not captured in the Landsat-derived surface water products, since the spatial resolution of Landsat products (30 m) is too coarse to capture those wetlands. This is typical of many surface water classifiers that are focused on deep open water, as they misclassify the highly variable spectral signatures of inundated wetlands [53]. Moreover, optical sensors struggled to capture wetlands covered by emergent vegetation. This study integrated Sentinel-1 SAR data into the high resolution (10 m) optical bands of Sentinel-2 to create a more robust classifier (Figure 9A). We also used a wider temporal window for the optical bands, which increased the number of observations over the study area. This allows our algorithm to minimize the effects of cloud covers, and identify the small wetlands by detecting frequently wet pixels. We performed an independent accuracy assessment on small and highly vegetated, and small non-vegetated wetlands. The results showed acceptable accuracy for both types of wetlands. We also compared the results with surface water maps derived from optical sensors (Table 6). Our algorithm performs better in identifying both large and small wetland water bodies than the Landsat-derived JRC and Sentinel-2-derived NDWI algorithms (Tables 4–6).

The wetland surface water was also evaluated in vegetated and non-vegetated wetlands. Visual observation shows that the small inundated wetlands contain more vegetation compared to larger and deeper water bodies. Comparing 2016 and 2017 wetland surface water maps reveals abrupt changes in emergent vegetation in small wetlands. These results agree with the findings of [54]. They reported that the small, ephemeral wetlands in the PPR experienced more vegetation change variability than larger, semi-permanent wetlands [54]. Large and deep water bodies can be easily detected by various remote sensing data. For example, [55] used Landsat time-series to create a global map of inland water dynamics. However, identifying small water bodies in the PPR is challenging due to the wetlands' size and strong potential for dense vegetation cover. This is very important, as the majority of wetlands in the PPR are small. This causes the surface water in potholes to be highly dynamic. The total surface water area calculated from the JRC product and our classification method was 294 km2 and 376 km2, respectively. Algorithms that miss surface water in these small wetlands will be biased, and misrepresent the hydrologic variability on the landscape. For example, small wetlands provide more foraging habitats for organisms that rely on shallow water.

Cloud computing and the advent of multisensor remote sensing data in the GEE have several advantages for large-scale and time-series analysis, such as monitoring wetlands dynamics [56]. The use of the GEE cloud computing platform is more convenient than traditional methods, considering its processing speed and ease of use [57]. As more machine learning algorithms and remote sensing data become available within the GEE platform, we expect remote sensing data processing to be simplified even further. Additionally, and unlike most supercomputing centers, GEE is also designed to help researchers quickly disseminate their results to other researchers and interested parties. Once an algorithm has been developed on the GEE, users can generate systematic data products or deploy interactive applications aided by the GEE's resources [25]. The fully automated workflow developed for this study allows us to refine the existing data and method, and rapidly

apply it to a broad geographical scale to generate estimates in new years. One of the disadvantages of using the GEE cloud computing platform is that it limits the number of field samples and input features. This is especially challenging when the analysis is applied to a large domain, which may reduce the efficiency of the implemented method.
