Monitoring and Analyzing the Seasonal Wetland Inundation Dynamics in the Everglades from 2002 to 2021 Using Google Earth Engine

: Inundation dynamics coupled with seasonal information is critical to study the wetland environment. Analyses based on remotely sensed data are the most effective means to monitor and investigate wetland inundation dynamics. For the ﬁrst time, this study deployed an automated thresholding method to quantify and compare the annual inundation characteristics in dry and wet seasons in the Everglades, using Landsat imagery in Google Earth Engine (GEE). This research presents the long-term time series maps from 2002 to 2021, with a comprehensive spatiotemporal depiction of inundation. In this paper, we bridged the research gap of space-time analysis for multi-season inundation dynamics, which is urgently needed for the Everglades wetland. Within a GIS-based framework, we integrated statistical models, such as Mann–Kendall and Sen’s Slope tests, to track the evolutionary trend of seasonal inundation dynamics. The spatiotemporal analyses highlight the signiﬁcant differences in wet and dry seasons through time and space. The stationary or permanent inundation is more likely to be distributed along the coastal regions (Gulf of Mexico and Florida Bay) of the Everglades, presenting a warning regarding their vulnerability to sea level rise.


Introduction
Inundation is the state of land when it is covered with water, either permanent or transient. Wetlands are the inundated or water-saturated lands characterized by various kind of vegetation [1], and inspecting inundation dynamics of wetlands is a critical first step towards understanding their hydroperiod, ecology, and functionality [2][3][4][5]. Wetland inundation dynamics is the key indicator of ecosystem functioning, and has a wide range of implications for water management and wetland governance [6][7][8]. However, it is often poorly understood owing to its complexity [7]. This is due to spatial and temporal heterogeneity, which creates uncertainties in the extent, timing, and duration of inundation. Considering the significant input of wetlands to local hydrology [9,10], wildlife habitat and flood mitigation [11], greenhouse gas emissions [12,13], and other environmental services, there is an urgent need for monitoring and analyzing the seasonal inundation dynamics through time and space.
There are several remote sensing techniques to quantify the inundation extents, such as spectral indices [14,15], regression trees [16,17], and the band ratio method [18]. However, these methods have been found to be less effective for mapping inundation pixels separately from non-inundation pixels. For example, using the indices technique, the normalized and modified normalized difference water index (NDWI and MNDWI, respectively), which were specially developed to extract waterbodies, have the least potential to classify inundation among the three methods [14]. Although the method relying on a regression tree is currently favored as it uses multi-scalar input data (such as thematic ancillary information), this model is still cumbersome and expensive as it requires a substantial The whole procedure to accomplish the materials and methodological part of this research is categorized into four broad sections, namely, study area selection (Section 2.1), data extraction and processing (Section 2.2), deploying the thresholding method to extract annual inundation characteristics from Landsat images (Section 2.3), and spatiotemporal analyses to investigate the spatial distribution and seasonal dynamics of inundation (Section 2.4).

Materials and Methods
A schematic framework was proposed in order to monitor and analyze the seasonal inundation dynamics of wetlands in the Everglades from 2002 to 2021 based on GEE (Figure 1). The whole procedure to accomplish the materials and methodological part of this research is categorized into four broad sections, namely, study area selection (Section 2.1), data extraction and processing (Section 2.2), deploying the thresholding method to extract annual inundation characteristics from Landsat images (Section 2.3), and spatiotemporal analyses to investigate the spatial distribution and seasonal dynamics of inundation (Section 2.4).

Study Area
Our study area lies in South Florida, United States (US), and comprises two unique and important national parks, the Everglades National Park and the Big Cypress National Park ( Figure 2). The Everglades region stretches from the southern rim of Lake Okeechobee (the second largest freshwater lake in the US) to Florida Bay [33]. This research

Study Area
Our study area lies in South Florida, United States (US), and comprises two unique and important national parks, the Everglades National Park and the Big Cypress National Park ( Figure 2). The Everglades region stretches from the southern rim of Lake Okeechobee (the second largest freshwater lake in the US) to Florida Bay [33]. This research selected the Everglades wetland as it is a very low-gradient environment [34], with an elevation range of approximately −0.65 to 5 m from the lowest to highest points [19]. The type of land cover encompassing the Everglades falls under the mosaics of open water and wet prairie with herbaceous plants, shrubs, and forested wetlands [35]. In addition, a major part of the Everglades is composed of wet prairie and sawgrass marsh [26]; this area has been named the "river of grass" for a long time [27]. According to National Land Cover Database [36] products 2019, the upper-most portion of the study area (Big Cypress National Park) is mostly dense vegetation (forest), whereas the wetlands is the most dominant land type covering the study area ( Figure 2). elevation range of approximately −0.65 to 5 m from the lowest to highest points [19]. The type of land cover encompassing the Everglades falls under the mosaics of open water and wet prairie with herbaceous plants, shrubs, and forested wetlands [35]. In addition, a major part of the Everglades is composed of wet prairie and sawgrass marsh [26]; this area has been named the "river of grass" for a long time [27]. According to National Land Cover Database [36] products 2019, the upper-most portion of the study area (Big Cypress National Park) is mostly dense vegetation (forest), whereas the wetlands is the most dominant land type covering the study area ( Figure 2).
With an annual mean temperature of 25 °C, the climate of South Florida is regarded as maritime subtropical, where summers are long and hot, and winters are mild and dry [37,38]. According to National Park Service, the Everglades has two seasons, namely, the wet season, which starts around the middle of May and continues through November, and the dry season, which ranges typically from December through April [39]. Hence, the wet season is characterized by intense rainfall and is longer than the dry season. The historic average rainfall of Florida from 1901 to 2000 was 136.27 cm [40].

Data Extraction and Processing
This research used remotely sensed data as a primary data source. Most of the satellite data processing was accomplished using GEE, unless otherwise stated. This research used images from three different Landsat sensors (Landsat 5 TM, Landsat 7 ETM+, Landsat 8 OLI), which were previously proven as providing effective source data for the study of inundation dynamics [28,41]. To retrieve the historical inundation characteristics, this study considered every tier 1 image, which were atmospherically corrected surface reflectance Landsat scenes captured from 2002-2021 in two seasons: the wet season ranges from May to November and the dry season from December to April. The Landsat sensors are capable of capturing scenes with 30 m spatial resolution, and they contain seven (Landsat With an annual mean temperature of 25 • C, the climate of South Florida is regarded as maritime subtropical, where summers are long and hot, and winters are mild and dry [37,38]. According to National Park Service, the Everglades has two seasons, namely, the wet season, which starts around the middle of May and continues through November, and the dry season, which ranges typically from December through April [39]. Hence, the wet season is characterized by intense rainfall and is longer than the dry season. The historic average rainfall of Florida from 1901 to 2000 was 136.27 cm [40].

Data Extraction and Processing
This research used remotely sensed data as a primary data source. Most of the satellite data processing was accomplished using GEE, unless otherwise stated. This research used images from three different Landsat sensors (Landsat 5 TM, Landsat 7 ETM+, Landsat 8 OLI), which were previously proven as providing effective source data for the study of inundation dynamics [28,41]. To retrieve the historical inundation characteristics, this study considered every tier 1 image, which were atmospherically corrected surface reflectance Landsat scenes captured from 2002-2021 in two seasons: the wet season ranges from May to November and the dry season from December to April. The Landsat sensors are capable of capturing scenes with 30 m spatial resolution, and they contain seven (Landsat 5), eight (Landsat 7), or eleven (Landsat 8) bands, including the SWIR band, which was used in this study. Pixels from each scene detected as a cloud or cloud shadow on the Landsat were masked using GEE's masking function. As a single sensor's images do not cover the examined time bound (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021), this research merged all image collections from Landsat 5, 7, and 8, and selected the SWIR band accordingly. Then, we used a gap-filling algorithm referred to as the Neighborhood Similar Pixel Interpolator (NSPI) method [42], which was developed based on the assumption that the same-class neighboring pixels have similar spectral characteristics when they are placed around the un-scanned pixels. With the help of the NSPI method, the masked pixels were filled using the median value for the pixel from the year before and after the scenes' date. The median value of the scenes for a particular year was used to create an annual composite for each year within the time limit.
To validate the accuracy of the classified annual inundation products, we used the continuous daily record of water depth maps from EDEN [31,43,44]. These datasets are preserved in 400 by 400 m grids, and all negative water depths were shifted to zero with the spatial interpolation method. To match the resolution with Landsat images and align with the time bound, we downscaled 400 m EDEN maps to a 30 m resolution, and derived the seasonal mean value (half yearly) from daily acquisition. Due to procedural difficulties and time constraints, we validated annual inundation products from the years of 2011, 2016, and 2021. Although EDEN maps are reliable local datasets for validating inundation products for our study area, they may have some anomalies due to the interpolation method they adapted during post-processing of the maps. Consequently, this research used a 10 cm threshold value of the water depth when we processed EDEN data for validation. The accuracy assessment of our classified inundation maps with EDEN water depth maps was based on pixel-to-pixel agreement of inundated and non-inundated areas. The total number of classified pixels of inundation and non-inundation was almost 4,300,000 in each participating year; however, we chose 7500 pixels using a random stratified sampling approach for validation. According to validation results, 2016 achieved the highest accuracy of 88%, followed by 84% in 2011, and 2021 had the lowest accuracy of 79%. The overall mean accuracy of annual inundation maps was 83.67% compared with EDEN maps, which is 12% greater than that of previous research [19].

Thresholding Method
We deployed the SWIR thresholding method, originally proposed by Wolski et.al [14], and further modified by Inman and Lyons [45], to transform the composite into annual inundation maps. Although we borrowed the synergetic version of the thresholding method from these two existing studies, we also modified and altered this model to fit in our spatial and temporal context. To apply this thresholding technique, we assessed and digitized inundated and non-inundated sample areas as polygons. When choosing typical polygons, we considered an area as inundated if it had contained permanent water (e.g., permanent swamp or main channels) based on Esri's high-resolution world imagery. According to the thresholding method, this study took the median value of the SWIR band for the inundated (SWIR inundated ) and non-inundated (SWIR non-inundated ) areas, and then it was applied for each annual image composite. The threshold (SWIR threshold ) for a specific composite was calculated using Equation (1): where SWIR inundated and SWIR non-inundated denote the median reflectance for the SWIR band in the inundated and non-inundated areas, respectively. SWIR threshold is the threshold output value, which defines inundated and non-inundated pixels within the study area. Pixel values below the SWIR threshold are classified as inundated, otherwise they are classified as non-inundated. As recommended by Wolski et al. [14], Equation (1) uses the multiplier of 0.3, which represents the required value for calculating the appropriate threshold to classify a pixel as inundated with a 50% inundation fraction. We also confirmed that using 0.3 as a multiplier in the threshold equation was appropriate due to different Landsat images used in this research.

Spatiotemporal Analysis
The spatiotemporal evaluation was conducted for seasonal inundation dynamics in terms of five aspects, namely, annual inundation distribution, overlapping and inundation dynamics, trend analysis, total inundation occurrence, and inundation variances. The detailed procedures and explanations are as follows: 1.
Annual inundation distribution: To obtain an overview of the inundation extent, it is essential to summarize the temporal change patterns of inundated areas. The areas were quantified using the pixelArea function in GEE and organized by year and season. To better compare and extract larger or smaller areas over the 20-year time span, this product is illustrated as bar charts for both dry and wet season.

2.
Overlapping and inundation dynamics: To inspect the temporal dynamics of inundation, an in-depth analysis was undertaken by overlapping and comparing the maps of inundation in the subsequent years [46,47]. A total of 19 comparison results were generated from the annual maps from 2002 to 2021, where inundated pixels in the previous year appearing as non-inundated pixels in the latter year were defined as a decrease event, and as an increase event in the reverse cases. Moreover, the pixels appearing as inundated/non-inundated in two consecutive years were regarded as an unchanged event. To better understand the temporal dynamics, the pixels were converted to area in km 2 , and summarized as bar charts with three events (decrease, increase, and unchanged) annually and seasonally. Beyond that, this dynamic comparison was performed between the annual maps of every four years (2002-2006, 2006-2010, 2010-2014, and 2014-2018), except 2018-2021, which was conducted with a 3-year interval due to data constraints. The result of this part (focusing on spatial distribution) was visualized on maps in three events (decrease, increase, and unchanged) for wet and dry seasons. The spatial patterns illustrated by the generated maps reveal significant spatiotemporal traits, in which pixel values assigned to 1 contributed to increase, −1 contributed to decrease, and 0 contributed to unchanged inundation.

3.
Trend analysis: An integrated nonparametric approach was used in this study to examine and assess if there is a trend within the inundation dynamics over time. This kind of approach is widely used in the GIS community [48][49][50][51][52][53][54][55][56]. Our research used the Mann-Kendal test [57,58] to investigate whether the seasonal inundation area follows any trend that is statistically significant (when p-value is less than 0.05). If it consists of any statistically significant trend, this study used Sen's slope estimator [59] to determine whether the trend runs through the downward or upward direction. To summarize the temporal trend of inundation dynamics, the annual inundation area from two seasons was applied, and the results were tabulated accordingly. 4.
Total inundation occurrences: An in-depth investigation was carried out based on the occurrence map by summing up the binary images from all of the 20 years. This was accomplished in GEE Code Editor using the sum function for image collection. As a result, the values of the pixels on the inundation occurrence map varied from 0 (without inundation occurrence) to 20 (always inundation occurrence). Based on the pixel value, the inundation occurrence can be categorized as very high (17)(18)(19)(20), high (13)(14)(15)(16), moderate (9-12), low (5-8), and very low (0-4). The spatial distributions of inundation can be obtained and compared in wet and dry season from this map. 5.
Inundation variance: As we had each year's inundation map from 2002-2021, we considered measuring the seasonal inundation dispersion from its mean inundation. The important insight of this variance map is the location of the hotspot of continuous inundation deviation within the study area. As with the total inundation occurrence map, this map was also produced in GEE Code Editor with the help of the variance function. The pixel values on the inundation variance map vary from 0 (without variance) to 0.25 (maximum variance). Based on the pixel value, the inundation variance can also be categorized as high (0.19-0.25), moderate (0.13-0.18), low (0.07-0.12), and very low (0-0.06). The inundation variance map was produced for season-to-season comparison.

Results and Discussion
The results and discussion are based on the analysis from the output maps and graphs, and they are elucidated and analyzed in this section. This section is divided into two broad sections: temporal dynamics of inundation (Section 3.1), and spatial dynamics of inundation (Section 3.2).

Temporal Dynamics of Inundation
Examining the annual inundation distribution over the years is a key analysis to retrieve the temporal change pattern of inundation. Figure 3 illustrates the extent of the inundated area for each year in wet and dry seasons from 2002 to 2021. The average inundated area per year is greater in the wet season (1484.81 km 2 ) than in the dry season (1373.33 km 2 ). However, analysis of the dispersion of inundation by season showed that the dry season signals more variability, with a standard deviation of 386.49 km 2 , compared to the wet season, which had a standard deviation of 292.88 km 2 . This comparison implies that the dry seasons from 2002-2021 were more active than wet seasons in terms of inundation dynamics. To further analyze annual inundation distribution, we explored the peak inundation years between the two seasons, and therefore present the longest ever time series of inundation extent in the study area. The observation from Figure 3a shows that the largest inundation extent (2076.31 km 2 ) in the wet season was recorded in 2015, followed by 1929.92 km 2 in 2002. The largest peak inundation was almost two times that of the smallest one (1006.69 km 2 ) recorded in 2012. The effect of the strongest El Nino episode in 2015 [60] could be one of the probable factors triggering this exceptionally large inundation extent compared with other years. However, other weaker El Nino events, such as the 2010 and 2020 episodes [61,62], presented a less remarkable change, which was also echoed in our study. The dry season (Figure 3b iance) to 0.25 (maximum variance). Based on the pixel value, the inundation variance can also be categorized as high (0.19-0.25), moderate (0.13-0.18), low (0.07-0.12), and very low (0-0.06). The inundation variance map was produced for season-to-season comparison.

Results and Discussion
The results and discussion are based on the analysis from the output maps and graphs, and they are elucidated and analyzed in this section. This section is divided into two broad sections: temporal dynamics of inundation (Section 3.1), and spatial dynamics of inundation (Section 3.2).

Temporal Dynamics of Inundation
Examining the annual inundation distribution over the years is a key analysis to retrieve the temporal change pattern of inundation. Figure 3 illustrates the extent of the inundated area for each year in wet and dry seasons from 2002 to 2021. The average inundated area per year is greater in the wet season (1484.81 km 2 ) than in the dry season (1373.33 km 2 ). However, analysis of the dispersion of inundation by season showed that the dry season signals more variability, with a standard deviation of 386.49 km 2 , compared to the wet season, which had a standard deviation of 292.88 km 2 . This comparison implies that the dry seasons from 2002-2021 were more active than wet seasons in terms of inundation dynamics. To further analyze annual inundation distribution, we explored the peak inundation years between the two seasons, and therefore present the longest ever time series of inundation extent in the study area. The observation from Figure 3a shows that the largest inundation extent (2076.31 km 2 ) in the wet season was recorded in 2015, followed by 1929.92 km 2 in 2002. The largest peak inundation was almost two times that of the smallest one (1006.69 km 2 ) recorded in 2012. The effect of the strongest El Nino episode in 2015 [60] could be one of the probable factors triggering this exceptionally large inundation extent compared with other years. However, other weaker El Nino events, such as the 2010 and 2020 episodes [61,62], presented a less remarkable change, which was also echoed in our study. The dry season (Figure 3b   To inspect the probable factors responsible for controlling the temporal distribution of inundation in the wet season, we looked at the historical dataset of rainfall occurrence and storm events during this time period (May-November) [63,64]. According to Figure 3a, the wet season in 2015 contained the highest inundated area; this is linked to the effects of Hurricane Danny and tropical storm Erka [65,66], which both occurred in August. According to the Florida Climate Center [67], the entire month of August 2015 witnessed high winds, hail, storm damage, and heavy rains across the state due to various convective and sea breeze thunderstorms. The long and substantial rainfall occurrence in South Florida is normally associated with tropical storms and hurricanes, which appeared between June and November [68]. Hence, the dry season in South Florida experiences less rainfall than that of the wet season unless there are exceptions during this time. According to Figure 3b, the dry season contains a less inundated area, as opposed to the wet season, in the overall observations. However, the dry season in 2007 witnessed increased inundation compared to the wet season in the same year, and this could be attributed to the off-season tropical storm Olga, which hit South Florida in December 2007 [69].
The results of overlapping comparison detailing the areas of three events, namely, increase, decrease, and unchanged inundation by year and season, are portrayed in Figure 4. Taking a closer look at the annual average area of decreased, unchanged, and increased inundation, it follows the ratios of 494:4514:532 in the wet season and 586:4327:627 in the dry season. Based on this ratio calculation, it is evident that the annual average area of inundation increased for the last 20 years in the dry season, and increased more than that of the wet season by approximately 18%. This is also true for the inundation decrease, where the annual average area in the dry season increased by 92 km 2 , which is roughly 19% greater than that of the wet season. The annual average area of unchanged inundation, in which inundation/non-inundation pixels appeared in two consecutive years, was found to be almost stable in both seasons, although the area of the wet season was 4% greater. The results of overlapping comparison mirror the findings from

p-Value
Sen's Slope Based on the annual inundated area (Figure 3), this study assessed the evolutionary trend of seasonal inundation dynamics using the Mann-Kendal test and Sen's slope estimator [48][49][50][51][52][53][54][55][56]. The results are presented in Table 1, including for the wet and dry seasons. In the Mann-Kendall test, the p-value less than 0.05 is sufficient to reject the null hypothesis (the trend does not exist) and accept the alternative hypothesis (the trend exists), whereas a positive value in Sen's slope estimator results in an upward trend and a negative value indicates a downward trend. According to Table 1, the inundation patterns of both wet and dry seasons fall under a statistically significant trend. Because Sen's slope estimator yielded a negative value in both seasons, the inundation trend follows a downward direction. Overall, the findings from trend analysis anticipate that, during the last 20 years (2002-2021), the whole study area witnessed a gradual decline in inundation, indicating less rainfall and precipitation occurrence in the Everglades and nearby regions.

Spatial Dynamics of Inundation
To extract the spatial scenario of inundation characteristics, the total inundation occurrence is a useful feature, and is illustrated in Figure 5. According to the map, the wet season (Figure 5a) presents more inundation features than the dry season (Figure 5b). Taking a closer look at the map, it is evident that a very high level of inundation clusters (pixel value [17][18][19][20] coincide along the coast of the Gulf of Mexico and Florida Bay. The area of the largest inundation cluster of the wet season (colored in blue) was approximately 0.39 km 2 , and the second largest area was roughly 0.33 km 2 . The dry season featured the largest cluster with an area of 0.35 km 2 , and the second largest area was roughly 0.29 km 2 . Overall, the largest inundation cluster in the wet season was greater than that of the dry season by roughly 10%, and all these large inundation clusters situated along the coast could be a scientific resource for future coastal research. Based on the season-to-season spatial dispersion of inundation features, the wet season contains larger inundation clusters than the dry season, indicating large drenched surfaces are prevalent from May to November. By contrast, high-level inundation clusters (pixel value [13][14][15][16] are witnessed mostly in the dry season, and these are predominantly distributed in the northern part of the study area. With the help of annual inundation maps from 2002 to 2021, the inundation variance map was generated and is demonstrated in Figure 6. The pixel values on the inundation variance map range from 0 (without variance) to 0.25 (maximum variance). Figure 6a displays the wet season, where high inundation variances (pixel value 0.19-0.25) are witnessed in the south-eastern part of the study area, which is adjacent to Miami-Dade County and the Water Conservation Area (WCA3B) [70]. Figure 6b, the dry season, also shows an almost similar distribution of inundation variances, with the exception of high-level variance in the upper region, and a slant cluster, located diagonally from the south-western part to the south-eastern part of the study area. One common attribute of inundation variances for both dry and wet seasons is that the two west coasts (Gulf of Mexico and Florida Bay) are characterized by less inundation variability, whereas the northern part, featuring Gator Lake and West Lake, and the south-eastern part, featuring Miami-Dade County and Tamiami canal, are identified as having a substantial change in inundation. With the help of annual inundation maps from 2002 to 2021, the inundation variance map was generated and is demonstrated in Figure 6. The pixel values on the inundation variance map range from 0 (without variance) to 0.25 (maximum variance). Figure 6a displays the wet season, where high inundation variances (pixel value 0.19-0.25) are witnessed in the south-eastern part of the study area, which is adjacent to Miami-Dade County and the Water Conservation Area (WCA3B) [70]. Figure 6b, the dry season, also shows an almost similar distribution of inundation variances, with the exception of highlevel variance in the upper region, and a slant cluster, located diagonally from the southwestern part to the south-eastern part of the study area. One common attribute of inundation variances for both dry and wet seasons is that the two west coasts (Gulf of Mexico and Florida Bay) are characterized by less inundation variability, whereas the northern part, featuring Gator Lake and West Lake, and the south-eastern part, featuring Miami-Dade County and Tamiami canal, are identified as having a substantial change in inundation.  With the help of annual inundation maps from 2002 to 2021, the inundation variance map was generated and is demonstrated in Figure 6. The pixel values on the inundation variance map range from 0 (without variance) to 0.25 (maximum variance). Figure 6a displays the wet season, where high inundation variances (pixel value 0.19-0.25) are witnessed in the south-eastern part of the study area, which is adjacent to Miami-Dade County and the Water Conservation Area (WCA3B) [70]. Figure 6b, the dry season, also shows an almost similar distribution of inundation variances, with the exception of highlevel variance in the upper region, and a slant cluster, located diagonally from the southwestern part to the south-eastern part of the study area. One common attribute of inundation variances for both dry and wet seasons is that the two west coasts (Gulf of Mexico and Florida Bay) are characterized by less inundation variability, whereas the northern part, featuring Gator Lake and West Lake, and the south-eastern part, featuring Miami-Dade County and Tamiami canal, are identified as having a substantial change in inundation. The spatial distribution of inundation information and its variability have a critical role to play in wetland management systems. As such, this research could be a potential resource to help the Comprehensive Everglades Restoration Plan (CERP), which aims to enhance freshwater storage, improve water quality and restore the natural water flow through the greater Everglades ecosystem [71]. By presenting a detailed depiction of total inundation occurrences ( Figure 5) and inundation variances ( Figure 6) maps, we attempted to enhance the spatial understanding and dynamic nature of inundation within the Everglades area. The results shown in Figures 5a and 6a indicate that a very high level of inundation occurrence (pixel level [17][18][19][20] and very low level of inundation variance (pixel level 0-0.06) almost overlap in the same coastal territory of the study area. This observation suggests that the coastal region of the Everglades stores permanent water from inland to a certain level throughout the wet season, therefore leaving less room for the dynamic nature of inundation. This kind of scenario indicates that the coastal area is vulnerable to the effects of sea level rise [43], which also means more tidal flats are converted to permanent water. Furthermore, the mangrove forest located along the southern Florida coast has been rapidly lost due to climate change [72], which calls for a meticulous awareness and further protection of the coastal environment of the Everglades.
This research investigated the possible factors that influenced the inundation occurrence and variance throughout the study area. To do so, firstly, we explored the elevation of the study area, which is regarded as one of the factors affecting the inundation events [73]. The long-term and recurrent inundation in both dry and wet seasons occurred at a low elevation of −0. 65-1.4 m [19,74]; this area is mostly located beside the coast of the Gulf of Mexico and Florida Bay. This accounts for roughly 30% of the total study area, where the results from Figure 5 show that areas belonging to low altitudes have a greater likelihood of being inundated throughout the year. Consequently, these areas need special attention of stakeholders and policymakers as flood and storm surges often strike in areas that have lower topographic elevation and downgraded streamflow [73]. The Big Cypress National Park, which constitutes the upper-most portion of the study area, has the highest elevation, ranging from 1.9 to 4.4 m [34], and therefore recorded less inundation compared to the lower part of the study area. Secondly, according to the LULC map of the study area (Figure 2), Big Cypress National Park is dominated by dense vegetation (forest), which accelerates the evapotranspiration process in this area, resulting in less inundation [74]. However, it is also possible that high forest cover can increase the precipitation rate over an area by importing evaporated moisture from the ocean [75].
To further analyze the spatial dynamic characteristics of inundation, we produced spatial distribution maps of overlapping comparisons between the inundation maps of every four years, and the comprehensive results are provided in Figures 7 and 8 for the wet season and dry season, respectively. In the wet season, the areas of unchanged inundation, which account for the largest portion of the study area, are 5170 km 2 in Figure 7c, 5251 km 2 in Figure 7d, and 5148 km 2 in Figure 7e. The majority of these large shares of unchanged inundation are detected in the central to upper parts of the study area, and this scenario indicates that the detected areas (colored in gray) have more stable or less active inundation activities than the other areas. In the dry season, Figure 8c-e also show the largest areas of unchanged inundation, which were 4812, 4965, and 4745 km 2 , respectively. However, in the dry season, the location of unchanged inundation clusters shifted towards the lower parts of the study area, and the upper parts were occupied by inundation increase (colored in blue) to a significant extent. Regarding inundation increase, Figure 7b shows the largest shared area of 863 km 2 , followed by Figure 7a, which covers an area of 773 km 2 in the wet season. In the dry season, Figure 8c shows the largest shared area of inundation increase (837 km 2 ), followed by Figure 8a (759 km 2 ). The largest shared area of inundation increase in the wet season is 3% greater than that of the dry season. Although the spatial distribution of inundation increase is mostly found in the central and northern parts of the study area (Figure 7a,b) in the wet season, the distribution of inundation increase in the dry season tends (Figure 8c,d) to be mostly in the southern part of the study area. When the area of decreased inundation is compared between dry and wet seasons (colored in red), the dry season experienced the largest inundation decrease (863 km 2 ), and it was distributed nearly throughout the study area (Figure 8b). The wet season experienced the highest inundation decrease, with an area of 539 km 2 , and this is also distributed roughly across the study area, with a higher concentration in the south-eastern part (Figure 7a). From the findings of this research, it is evident that the northern portion of the study area, which features Big Cypress National Park, is dominated by a dynamic increase in inundation. That means that this area had an active participation with inundation dynamics throughout the year. Another active area is the south-eastern portion of the study area, where it showed a significant inundation decease in the last 20 years, which is a warning sign of an environmental crisis.
Although it is hard to pinpoint the exact causes of inundation differences, the existing literature [76] suggests the differences may arise from short-term agents and long-term consequences. Short-term agents causing coastal inundation include storm surges and hightide flooding, while long-term causes result in relative changes in local and/or global sea level. As one of the driving forces of coastal flooding [77], sea level rise inundates low-lying wetlands and dry lands, causes shoreline erosion, and contributes to the influx of saline water into estuaries and nearby groundwater aquifers. Moreover, coastal infrastructures are more vulnerable to the damage arising from higher sea level rise caused by storms. Based on several previous studies, including those of the US Environmental Protection Agency (EPA) [78][79][80][81][82], inundations have become frequent along the US coastline since the 1950s, and the resultant sea level rise is affecting human activities in coastal areas. In the past decade (2011-2020), the coast of Florida facing the Gulf of Mexico experienced extended inundation that was above the local threshold limit for flooding along US coasts. In particular, the two sites of St. Petersburg and Cedar Key, within or adjacent to our study area, were inundated for an average of 15 days per year due to coastal flooding [80].

Conclusions
This research deployed an efficient thresholding of the SWIR band method and produced annual inundation maps from long-term time series of Landsat images (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021) in the Everglades. Our study conducted several spatiotemporal analyses on the produced inundation maps to investigate the dynamic inundation in terms of annual and seasonal aspects. The inundation results highlight the significant differences in wet and dry seasons in time and space, where the wet season was found to be more stationary and less active than the dry season in terms of inundation dynamics. Although the permanent inundations are more distributed along the coastal regions (Gulf of Mexico and Florida Bay) of the Everglades, the results highlight their vulnerability to an environmental crisis, and therefore call for a higher level of precaution and protection.
The detailed spatiotemporal background and season-to-season comparisons with analytical reasoning of inundation dynamics are the key contributions of this study. In this paper, we fulfilled the research gap of space-time analysis for multi-season inundation dynamics, which was urgently needed for the Everglades wetland. We anticipate that these periodic and accurate inundation products, along with comprehensive spatiotemporal depictions, could be used by land managers, researchers, and other stakeholders who require access to high-resolution inundation maps. With this research, we showed that GIS-based quantitative analyses are more powerful than conventional techniques, and therefore could be used to better inspect and accurately portray the dynamic changes of geographic events on large spatiotemporal scales. Although this study implemented the thresholding method to depict the seasonal inundation dynamics in the Everglades, this could be applied with minimal adaptation in other wetland studies. Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://developers.google.com/earth-engine/datasets/ (accessed on 13 November 2022).