Comprehensive Landsat-Based Analysis of Long-Term Surface Water Dynamics over Wetlands and Waterbodies in North America

Abstract

Wetlands are considered one of the most valuable ecosystems around the world and provide numerous environmental services, including water purification, flood protection, and habitat for a variety of species. Wetlands loss is an increasing trend due to anthropogenic activities and natural processes. As such, spatial knowledge regarding the extent and dynamics of surface water is demanding for wetland conservation and protection. The Landsat program, with five decades of historical Earth observation data, has a unique advantage for monitoring wetland surface water changes and dynamics with 30 m spatial resolution. We monitored 266 Ramsar wetland sites in North America for the past 40 years using the open-access Landsat data within the Google Earth Engine cloud computing platform. Landsat Collection 2 Level-2 surface reflectance products were preprocessed and cloud-screened, and a time series of spectral bands and indices were created. The unsupervised Dynamic Surface Water Extent method classified each image into water classes with different confidence levels. An average overall agreement of 92,97% and an average F-score of 96.31% were achieved in this study. Water occurrence maps, in addition to inundation class and change maps, were created for the entire North America, and quantified spatial information was calculated for Ramsar wetland sites.

RÉSUMÉ

Les zones humides sont considérées comme l’un des écosystèmes les plus précieux au monde et fournissent de nombreux services environnementaux, notamment la purification de l’eau, la protection contre les inondations et l’habitat d’une variété d’espèces. La perte de zones humides est une tendance à la hausse en raison des activités anthropiques et des processus naturels. En tant que tel, les connaissances spatiales concernant l’étendue et la dynamique des eaux de surface sont exigées pour la conservation et la protection des zones humides. Le programme Landsat, qui s’appuie sur cinq décennies de données historiques d’observation de la Terre, présente un avantage unique pour surveiller les changements et la dynamique des eaux de surface des zones humides avec une résolution spatiale de 30 m. Nous avons suivi l’évolution au cours des 40 dernières années de 266 sites Ramsar de zones humides en Amérique du Nord à l’aide des données Landsat en libre accès et de la plate-forme Google Earth Engine. Les produits de réflectance de surface de niveau 2 de la collection Landsat 2 ont été prétraités et les nuages identifiés, et des séries chronologiques de bandes spectrales et d'indices ont été créées. La méthode non supervisée de l’étendue dynamique des eaux de surface a classé chaque image en classes d’eau avec différents niveaux de confiance. Une concordance globale moyenne de 92.97% et un score F moyen de 96.31% ont été obtenus dans cette étude. Des cartes d’occurrence de l’eau, en plus des cartes des classes d’inondation et des cartes de changement, ont été créées pour l’ensemble de l’Amérique du Nord, et des informations spatiales quantifiées ont été calculées pour les sites Ramsar de zones humides.

Introduction

Various critical and important ecosystem services are provided by wetland ecosystems, including erosion reduction, carbon storage, flood protection, groundwater recharge, wildlife habitat, and water filtering and purification (Mahdianpari et al. 2021; Zhao et al. 2019). Despite these benefits, wetland environments have been altered or destroyed globally due to anthropogenic reasons, including industrial development, dams, urbanization, human disturbance, climate change, and agriculture (Davidson and Davidson 2014; Hemati et al. 2023). Retraction of floodwaters across the floodplain and the predictable advance maintained the floodplain wetlands (Junk et al. 1989). Herbaceous plants and grasslands tend to dominate the wetland areas, and poorly drained mineral soils underlain these ecosystems, causing seasonal or permeant water inundation. Rainfall and snowmelt events can change the extent and level of floodplain water levels rapidly and frequently (Lyon et al. 2010), which initiate a series of biological processes (Allen et al. 2020).

Wetland areas are defined by the United States Army Corps of Engineers (Environmental Laboratory 1987) as “those areas that are inundated or saturated by surface or groundwater at a frequency and duration sufficient to support, and that under normal circumstances do support, a prevalence of vegetation typically adapted for life in saturated soil conditions.” Using this definition, the small and large waterbodies are included in this study (Halabisky et al. 2016). Surface water is one of the most important parts of wetland ecosystems, and information regarding water inundation and changes is vital for wetland conservation. Hence, wetland surface water dynamic has been the subject of various studies investigating the dynamics of wetlands through monitoring the surface water (Halabisky et al. 2016; Jin et al. 2017; Wu et al. 2019). Spatial and temporal knowledge of the Earth’s surface waters is crucially important for different levels of policy decisions, from local to the national or global level, on the quality of life, wildlife, and the environment (Hemati et al. 2022). From areas that are permanently covered by water to non-water dry land, the location and frequency of surface water can vary. This region with temporary water inundated areas between those extreme permanent water and non-water areas experiences floods periodically (Olthof and Rainville 2022).

Precipitation and temperature trends have been altered with climate change across many regions, causing changes in inundation frequency and extent (Owen et al. 2020). Evapotranspiration and precipitation patterns are expected to alter in the northern hemisphere due to climate change associated with warmer air temperatures (Feulner et al. 2013; Kouhgardi et al. 2022). Wetlands hydrology may be affected by such changes by reducing available ground and surface water sources (Rokaya et al. 2020), and this may further stress water quality and ecosystem biota (Leach et al. 2020). Therefore, there is an increasing demand for knowledge and research into surface water and wetland monitoring and mapping to expand our understanding of long-term wetland changes in the past decade and their vulnerability. Adopted on 2 February 1971, the Ramsar Convention is a global intergovernmental agreement between multiple nations (i.e., 170 countries as of now) with the purpose of conservation, sustainable use, and restoration of wetlands. This multilateral environmental agreement covers different aspects of wetlands, including detailed reports about wetland types, vegetation, waterbirds, fishes, and ecosystems. It also consists of management plans for the wise use of natural resources (Ramsar Convention 2016). However, the Ramsar Convention lacks sufficient spatial information, missing wetland extent, maps, and water dynamics (Ramsar Convention 2016).

Considering the coverage and temporal resolution of different Earth Observation (EO) data, remote sensing is one of the most cost-effective and accurate tools to monitor the dynamics of wetlands and surface waters (Hopkinson et al. 2020). The Landsat mission has provided continuous science-grade open-access satellite images for the past 50 years with global coverage and radiometric consistency in time and space (Hemati et al. 2021). This data has paved the way for policy support activities and different environmental applications, including wetland mapping and monitoring (Wulder et al. 2022). Other remote sensing instruments, such as the Sentinel satellite constellation, provide more frequent optical and radar data with higher spatial resolution. However, they lack the historical archive and spectral consistency of the Landsat program (Wulder et al. 2019). Only in recent years has applying remote sensing applications on a national or continental scale become feasible with almost no cost due to the availability of geospatial cloud computing platforms, such as Google Earth Engine (Gorelick et al. 2017). The GEE platform overcomes large-scale data handling challenges and provides the resources required to store and process a massive amount of EO data.

Dynamic maps have been produced from a time series of water masks based on the EO data archive to depict and monitor the frequency, location, and changes of surface water. From regional to continental and global scales, maps have been produced in different spatial resolutions over multiple decades. Dynamic surface water maps usually depict the percentage or probability of water present in a time series of valid (not cloudy) observations for a pixel as water occurrence or inundation frequency (Mueller et al. 2016). However, knowledge regarding the time of the presence of surface water and inundation occurrence is not clearly presented by water occurrence maps. For example, a pixel with 50% water occurrence can be from an area where water is observed only half of the year due to seasonality (transition between dry and wet states intra-annually) or from an area where water is detected on one-half of the period and that area dried on the other half (permanent change; Olthof and Rainville 2022). To describe the nature of dynamic surface water and provide complementary information in addition to water occurrence maps, studies provide change maps (Pekel et al. 2016), post-classification comparison (Prigent et al. 2020), and logistic regression inundation frequency trend analysis (Olthof and Rainville 2022). Decent reference data from available features on the Earth (ground truth) and knowledge of associated spectral properties are required as the training data for supervised classification or as the validation data for the accuracy assessment of the produced maps. Considering logistical challenges and the cost of collecting reference data, especially for large-scale studies, ground truth data may not be available, or the amount of it is sub-optimal for accurate results. Additionally, surface water, as a dynamic target, does not have a consistent spectral response in optical remote sensing due to imaging conditions, depth, ice cover, chlorophyll, and turbidity (Dekker et al. 1997).

In one of the most significant and notable studies, over four million Landsat scenes have been used to map and monitor global inland and coastal surface water from 1984 to 2015, with the benefit of ongoing yearly updates until 2021 (Pekel et al. 2016). This impressive data includes occurrence, change magnitude, seasonality, recurrence, max extend, and a categorical classification of transitions, and it is publicly available on GEE and also on the Global Surface Water Explorer online portal provided by the European Commission (European Commission 2020). Other studies also used Landsat data for Global inland water dynamics (Pickens et al. 2020), seasonal dynamics (Pickens et al. 2022), and regional applications (Halabisky et al. 2016; Jin et al. 2017). Besides Landsat, MODIS (Alonso et al. 2020), Sentinel-2 (Heidi van et al. 2022; Yang et al. 2020), Sentinel-1 (Tang et al. 2022), and RADARSAT-2 (Zhao et al. 2014) were used to monitor dynamics of surface water, but they lack the temporal or spatial coverage of Landsat data, or they had a coarse spatial resolution.

Previous surface water dynamic studies were more focused on water bodies; hence, their algorithm was developed to detect water and non-water areas (Olthof and Rainville 2022; Pekel et al. 2016). However, in complex wetland areas, there is a wet region between open water areas and upland areas that is challenging to describe with binary water maps. Hence, in this study, a multiclass water detection method was proposed for dynamic surface water mapping. In this way, there is more flexibility to map and monitor water dynamics, especially on waterbodies in the wetland regions, with water classes with different confidence levels. Studies on wetland surface water dynamics previously were more focused on waterbody area (Halabisky et al. 2016) and maximum extent (Wu et al. 2019). Moreover, the long-term wetland water dynamic studies were limited to a single wetland site or small region (Hopkinson et al. 2020; Jin et al. 2017) or had limited temporal coverage (Amani et al. 2022; Battaglia et al. 2021). Therefore, a modified version of the Dynamic Surface Water Extend method, which was limited to the U.S. before, was adopted in the GEE. The modification was done to provide a scalable method to monitor surface water dynamics, even in higher latitudes, and to analyze time series, which were not implemented and applicable before. Additionally, the launch of Landsat-9 and reprocessing of Collection 2 data provides a new opportunity to produce dynamic surface water maps of wetlands with improved geometric, radiometric, and atmospheric correction of Landsat time series than previous studies.

The main focus of this study is on Ramsar wetland sites. However, the smaller wetland sites, lakes, rivers, and various forms of surface water were mapped in this study to provide spatiotemporal information regarding different types of waterbodies in addition to wetland surface water dynamics. The large-scale long-term maps produced in this study highlight spatial and temporal patterns in wetland surface water changes that are unique and provide more advanced information for local and governmental decision-makers. Briefly, the objectives of this study are:

• Monitoring surface water dynamics for wetland areas and waterbodies, using 40 years of Landsat Collection 2 data at a large scale over North America.

• Implementing a modified version of the Dynamic Surface Water Extent algorithm in the GEE to use Collection 2 time series and producing water dynamic maps in different latitudes.

• Change detection and classification of surface water dynamics and calculating quantitative reports at a wetland site level for registered wetlands at Ramsar Convention.

Therefore, in this study, the historical archive of Landsat and GEE cloud computing were leveraged to produce dynamic water maps of North America with temporal and spatial information regarding surface water extent, inundation frequency, and changes in wetland areas and waterbodies.

Materials and methods

We used historic EO data from the Landsat program, made publicly available by the United States Geological Survey (USGS). This data was accessed within the Google Earth Engine (GEE) cloud computation platform. An overview of the flowchart of this study has been illustrated in Figure 1.

Figure 1. Flowchart of the method used in this study for wetland’s surface water monitoring.

Study area and wetland sites

North America is a continent on the North American Tectonic Plate, bordering the Atlantic Ocean in the east, the Arctic Ocean in the north, and the Pacific Ocean in the west. Canada, the United States, and Mexico are the biggest countries in North America and, therefore, the main focus of this study. However, some other smaller countries such as Nicaragua, Honduras, Cuba, Guatemala, Dominican Republic, Haiti, Belize, Bahamas, Jamaica, Dominica, St. Lucia, Antigua and Barbuda, Barbados, and several Caribbean islands, create our study area, totaling 34 countries and territories. The study area covers around 24,939,249 square kilometers, about 16.5% of Earth’s land area, and is home to more than 592 million people (United Nations 2019). To manage the big data processing, the study area was divided into 15° × 15° grids (Figure 2), creating a network of 32 smaller parts, covering North American countries between 10°–85°N and 65°–170°W.

Figure 2. The grid network covers the study area and Ramsar Sites.

This area contains 220 Ramsar wetland sites covering about 23,628,375 hectares. Mexico has the most Ramsar sites in this region, with 142 sites, followed by the United States and Canada, with 41 and 37 sites, respectively. There are also 46 Ramsar Sites located in the study area region from other smaller countries, increasing the number of wetland sites to 266, wherein 137 of them are included in Ramsar Criterion 1, meaning that they are considered internationally important, as they contain a natural or near-natural wetland type that is unique, representative, or rare in the appropriate biogeographic region.

Earth observation data processing

Landsat program is the longest-running EO program, imaging the Earth from satellites for the past five decades, since 1972. This unique historical open-access archive contains an unbroken series of multispectral science-grade observations with medium spatial resolution. The number of Landsat scenes used from different Landsat platforms is illustrated in Figure 3.

Figure 3. The number of Landsat images from each Landsat platform.

Landsat-4, launched in 1982, was the beginning of the new generation of Landsat satellites with new technological developments, followed by Landsat-5, which holds the record for longest working satellites (29 years) and provides 397,601 scenes in this study. It should be noted that data from Landsat-1, -2, and -3 were not used in this study due to different spatial resolutions, lack of SWIR bands, and different scene reference systems. In June 2003, a Scan Line Corrector (SLC) failure occurred for Landsat-7, and since then, the data has been delivered with gaps. However, due to pixel-based analysis in this study, there is no issue using this data, and gaps are considered as other masked parts, such as clouds. On 27 September 2021, the continuity of the Landsat program was secured with the launch of Landsat-9 with improved imaging quality. The number of Landsat scenes used in this study for each year is depicted in Figure 4.

Figure 4. The number of Landsat images accessed for each year from 1982 to 2022.

In 2013, the highest number of images were available, with 40,456 scenes. Between 1984 and 2013, 13,000 to 17,000 scenes were available each year, with a peak of over 30,000 scenes between 2000 and 2002. At this time, Landsat-5 perfectly covered the failure of Landsat-6. With the launch of Landsat-8 in 2013, almost 30,000 scenes were made available each year due to increased onboard storage and transfer rates (i.e., an average of 740 scenes per day globally). The number of Landsat scenes available for each month is shown in Figure 5.

Figure 5. The number of Landsat scenes for each month of the year.

For capturing and monitoring the surface water dynamics of wetlands through the dry and wet seasons, a minimum amount of data for each month of the year is required. Considering the 16-day temporal resolution of Landsat satellites at the equator, the monthly number of images was assumed to be consistent throughout the year. However, due to the occurrence of a big part of the study area in higher latitudes (Canada and Alaska), a seasonal pattern was discovered for the number of available Landsat images each month in the GEE. The highest number of Landsat scenes was available for July, with 105,147 scenes, and the lowest was available in December, with 30,625 scenes. Considering the increased amount of data in some months of the year, the number of cloud-free scenes (i.e., with less than 20% cloud coverage) did not change with the same trend due to more cloud coverage in those months. The spatial distribution of cloud-free Landsat data over North America for the past four decades is illustrated in Figure 6.

Figure 6. Spatial distribution of Landsat data with cloud coverage less than 20% over North America between 1982 and 2022.

In this study, we used about a million (i.e., 928,748) Landsat scenes available from 1982 to 2022 in North America. As recommended by previous studies (Olthof and Rainville 2022; Tulbure and Broich 2013) and based on the trial-and-error preliminary results, highly cloudy images were excluded. Based on monthly analysis (Figure 5), we found enough observations to monitor the dynamics of surface water, only using less cloudy images. Excluding cloudy images not only increased the quality of the time series but also required less computational time. The second reprocessing of the Landsat archive, known as Collection 2, was produced in September 2020 by the United States Geological Survey (USGS) to achieve the highest consistency within the Landsat archive. Level-2 Surface reflectance data were used in this study to ensure radiometric consistency within the time series to minimize the errors related to different sensors and atmospheric conditions. Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm (version 3.4.0) for Landsat-4, -5, and -7 and Land Surface Reflectance Code (LaSRC) for Landsat-8 and -9 were used to produce surface reflectance data. The best achievable geolocation accuracy was provided for time-series by L1TP (precision and terrain corrected) standard with a geolocation root-mean-squared error (RMSE) less than 12 meters for Tier 1 data, which was used in this study. These images were accessed from Landsat-4, -5, -7, -8, and -9, available within GEE, and scenes with cloud coverage of less than 20% were filtered. For each image, we selected the visible RGB, Near-Infrared (NIR), Short-wave Infrared (SWIR1,2), and Pixel quality (QA_PIXEL) attributes, and time series of images were constructed. Bitmask produced by the CFMASK algorithm (Foga et al. 2017) was used to create a new band named "clouds" and mask clouds and dilated clouds (Bit 1 and 3), cloud shadows (Bit 4), snow (Bit 5), and cirrus clouds (Bit 2, available in Landsat-8 and -9).

In addition to Landsat data, the ALOS World 3D (AW3D30) Digital Surface Model (DSM) was used to include the geometric features. This dataset was preferred over the SRTM DEM that was mostly used in other studies and original Dynamic Surface Water Extent products by USGS since the SRTM does not cover most parts of Canada and Alaska (i.e., above 60°N) in the GEE. AW3D DSM was produced by stereo pairs of optical images and an image matching process, which has a 30 m spatial resolution (1 arcsec mesh) that is compatible with Landsat spatial resolution.

Dynamic Surface Water Extent

Dynamic Surface Water Extent (DSWE) was originally developed by USGS (Jones 2019) for Landsat Level-3 science products and is only available for conterminous U.S., Alaska, and Hawaii. This algorithm is based on several spectral indices that have been summarized in Table 1.

Table 1. Spectral indices have been used for the DSWE algorithm.

Spectral Index	Full Name	Formula
NDVI	Normalized Difference Vegetation Index	(NIR - Red) / (NIR + Red)
MNDWI	Modified Normalized Difference Wetness Index	(Green - SWIR1) / (Green + SWIR1)
MBRSV	Multi-band Spectral Relationship Visible	Green + Red
MBSRN	Multi-band Spectral Relationship Near-Infrared	NIR + SWIR1
AWEsh	Automated Water Extent Shadow	Blue + (2.5 * Green) + (−1.5 * MBSRN) + (−0.25 * SWIR2)

DSWE indices are designed to be sensitive to water regardless of the spectral reflectance change of water due to depth, aquatic vegetation, bottom soil characteristics, and its constituents. Using the sun azimuth and elevation, available in the metadata of Landsat scenes, and ALOS DSM elevation, slope, and hillshade were calculated for each image. The DSWE algorithm is based on five spectral tests applied to spectral bands and indices. These tests are designed for the unsupervised classification of each image into four different water-inundated classes. These tests are rigorously adjusted for Landsat data and require spectral reflectance data as a physical value consistent from time to time and place to place for the same land cover. Landsat Collection 2 Level-2 surface reflectance data contains scale factor and offset value for each band that was first applied to digital number values, and then the thresholds and spectral tests were applied to images. These thresholds are defined by a linear spectral mixture model, and the default values suggested by the USGS are used in this study.

There have been several studies examining the compatibility between reflectance and spectral indices obtained from the Enhanced Thematic Mapper (ETM) and the Operational Land Imager (OLI) (Flood 2014; Roy et al. 2016). At-surface reflectance calibration found 3–5% discrepancies between the two sensors, depending on the calibration method (Mishra et al. 2014). Notable disparities were identified between OLI- and ETM-derived values of the Normalized Difference Vegetation Index (NDVI) and Modified Normalized Difference Water Index (MNDWI) at low and high values, respectively (Li et al. 2013). Nevertheless, the DSWE model utilizes rule thresholds at elevated values for NDVI and extremely low values for MNDWI, rendering it impervious to the observed discrepancies. Inferences made from vegetation indices using in-situ data resampled to OLI/ETM spectra (She et al. 2015) suggest that surface reflectance spectra derived from ETM and OLI exhibit adequate correlation for the purpose of longitudinal studies (Jones 2019).

To apply spectral tests, a 5-bit bitmask was created, in which each bit represents the result of each test. For the first test on the Modified Normalized Difference Wetness Index (MNDWI), we masked the area with a value above 0.124. The second test masked the area where the Multi-band Spectral Relationship Visible (MBSRV) is more than the same for near-infrared (MBSRN). The third test selects the area where the Automated Water Extent Shadow is greater than zero. The fourth and fifth tests are specified for partial surface water (PSW) and have multiple conditions. PSW1 (Test-4) highlights the area where MNDWI is greater than −0.44, SWIR1 is less than 900, NIR is less than 1,500, and NDVI is less than 0.7. The last test (PSW2) includes five conditions, including MNDWI higher than −0.5, SWIR1 less than 3,000, SWIR2 less than 1,000, NIR less than 2,500, and Blue less than 1,000. Each pixel is then classified into non-water or four water classes: open water (high confidence), open surface water (moderate confidence), partial surface water, and water or wetland (low confidence). Considering the general low information in a shadowed area, the hillshade was used to mask areas with values less than 110. Slope values are also used to mask areas with high slopes (16.7° for the first two open water classes) to reduce commission error.

Water occurrence and change detection

To collect information about water inundation frequency, the spatial coverage of surface water, and the dynamics of wetlands, water occurrence maps were produced using the steps shown in Figure 7.

Figure 7. The flowchart for producing water occurrence and change detection maps.

After the construction of the time series for each period, each image was preprocessed, meaning cloud, cloud shadow, and snow were masked, and spectral indices and geometric features were extracted. Each scene is then classified using the DSWE unsupervised method. By selecting high and moderate confidence level open water classes, a time series of water and non-water was created for each pixel based on clear observations. The number of water-inundated pixels is normalized based on the number of clear observations for each pixel, and water occurrence maps are created. To add complementary information regarding the dynamics of surface waters and wetlands and to identify long-term changes, change maps were produced based on differential change detection methods and water occurrence maps of two different periods. Based on a threshold on the water occurrence map, the inundation classes have been classified into permanent and temporary water. In addition, to monitor long-term changes in Ramsar wetland sites and surface water dynamics, water loss, and water gain maps are produced using change magnitudes. Statistical analysis based on inundation class and change categories is also conducted for each wetland site in this study. Finally, to evaluate the results, the European Commission’s Joint Research Center (JRC) Yearly Water Classification History v1.4 (Pekel et al. 2016) has been used in this study. This data has been made accessible with the partnership of Google and UN Environment, is freely available in the GEE data catalog, and has been updated and reproduced until 2021. Annual water masks from JRC were compared to yearly water masks created by DSWE. Overall agreement and F-score were calculated for each year to show the performance of the proposed method throughout the time series.

Results

The feasibility of the Landsat historical archive for monitoring the long-term dynamics of surface water and changes in wetlands and waterbodies using a cloud computing platform was shown in this study. Using 32 gridded sub-regions, the 25 million square kilometers of North America were mapped with 30 m spatial resolution. We applied a 0.3° buffer to official border lines of North American countries to include information about coastal wetlands and shoreline changes in this study. Figure 8 shows the water occurrence map of North America from 1982 to 2022, which was produced in this study.

Figure 8. The water occurrence map of North America from 1982 to 2022 with 30 m spatial resolution.

All types of surface water were mapped in this study, ranging from shallow open-water wetlands to rivers and large lakes. Lake Erie, Huron, Superior, Ontario, and Michigan, known as the Great Lakes of North America, are the most extensive bodies of water in the water occurrence map, located on the border of Canada and the United States. These sea-like freshwater lakes are interconnected by the Great Lakes Waterway, which enables travel and shipping between them. Lake Superior is the biggest body of water detected in the results, with an estimated area of over 82 thousand square kilometers. Located in the Northwest Territories of Canada and in the boreal forest, Great Bear Lake is the biggest body of water mapped in this study, which is entirely in Canada. Lake Winnipeg, Grate Slave Lake, and Lake Athabasca are the other big water bodies in Canada. In addition to water occurrence maps and change results of the Quill Lakes, the summer composite of Landsat satellites for every five years, along with corresponding unsupervised classification results, have been demonstrated in Figure 9 to visually assess the performance of the DSWE method in time series. Quill Lake was chosen based on the intense changes detected by previous studies (Olthof and Rainville 2022) as an excellent example to show the capabilities of the proposed method on a fine scale.

Figure 9. Results of DSWE, water occurrence, and change map for the Quill Lakes, Saskatchewan, Canada (51.90 N, 104.36 W).

This wetland complex was designated an internationally important wetland via the Ramsar Convention in 1987. Initially, this wetland complex consisted of Big Quill Lake, Middle Quill Lake, and Little Quill Lake. Three separated water bodies can be seen in the early images. The visual comparison of true-color images and DSWE unsupervised classification results show that this method was perfectly able to separate open water and surrounding partial waters. As seen, there was a gradual change in the water level from 1985, and the water level has increased significantly since 2010. In 2015 and 2020, all three lakes were connected due to the permanent flooding of this wetland complex in the summer. The water occurrence map from 1982 to 2002 illustrates the extent of permanent and temporary inundated areas. Following the trend, the water occurrence map between 2002 and 2022 depicted the surface water dynamics of the wetland complex. To capture this long-term change in the wetland water occurrence, we produced the change map, which shows the magnitude of the changes (i.e., −1 for water loss, +1 for water gained).

To illustrate the novelty and improvements of surface water detection in this study, the water classification and occurrence maps produced in this study are compared to the products of the JRC (Pekel et al. 2016). Three examples in different latitudes (12, 29, and 52 N degrees) were chosen to show the performance of the proposed method in different environmental conditions. A decision tree is used by the JRC to classify each image into water and non-water areas, and yearly classification maps are created based on the water occurrence maps and temporal information. However, in the proposed approach, each image is classified into four water classes, and water occurrence maps can be produced with higher flexibility. In the JRC classification map, wetland surface water is classified under seasonal or permanent waters based on the frequency of the water inundation and temporal data. In the DSWE maps, surface water is classified into four water classes based on the confidence level. In this way, regardless of the temporal information, surface water classification can be done for a single Landsat scene or monthly, seasonal, and yearly composites.

Mitchell Lake and Laguna Monte Galan show two lakes surrounded by wetland vegetation in the U.S. and Nicaragua, respectively. Despite the permanent inundation of water, JRC water classification falsely classified central parts as seasonal water and failed to classify boundary areas where wetlands occur. The JRC water occurrence in these two lakes also ranges from 30 to 40, indicating the failure of the algorithm used in the detection of water in a great part of the time series. However, in the DSWE classification map, water bodies and boundary areas are perfectly classified with different water classes. The advantage of water detection with multiple classes is illustrated, and the flexibility of the proposed method maps the inundation frequency accurately. In a more severe example, Sullivan Lake from Canada is picked, where the JRC classification map is very noisy and completely fails to map the lake. The failure of the JRC method to identify water inundation is more obvious in the water occurrence map, where the values are around 10%. This clearly indicated that this permanent water body was misclassified in most of the time series. The flexibility and novelty of the proposed method are again shown in the DSWE classification map, where the water is detected with moderate confidence. This multiclass system and water detection with different confidence levels overcame the issue in the JRC product, and the water occurrence map was produced accurately. Despite the advantages of the produced products in this study over JRC products, the JCR dataset was used to measure the overall performance of the proposed method since it covers North America and is available in the GEE for the desired period.

The validation results have been illustrated in Figure 11. The overall agreement was calculated between the produced annual maps by DSWE in this study and the JRC yearly water classification product to show the performance and quality of the results throughout the time series. Despite the advantage of the produced studies (Figure 10), the JRC product proved to be a golden standard and has been used as a reference dataset in previous surface water dynamic studies for evaluation or comparison (Olthof and Rainville 2022; Wu et al. 2019; Zou et al. 2021). The overall agreement shows the number of samples detected as water in both products from all samples. The consistent performance of the proposed method was shown, and an average overall agreement of 92.97% was achieved. Furthermore, the harmonic mean of recall and precision values from the produced maps were calculated as the F-score. This validation index shows the producer’s and user’s accuracy of the water detection in a single measurement. The average F-score was 96.31%, showing acceptable water class detection from 1984 to 2021.

Figure 10. Visual comparison of the Landsat annual median composite image of 2020, water classification, and water occurrence maps from JRC and this study. (Mitchell Lake, San Antonio, Texas, USA; Laguna Monte Galan, León, Nicaragua; and Sullivan Lake, Alberta, Canada).

Figure 11. Overall agreement and F-score of the produced results in different years obtained from JRC yearly water classification product.

General information about wetland type, Ramsar Criteria, and management plans is available from the Ramsar Sites Information Service online portal. However, this information is generally outdated, only gives an overview of each wetland site in North America, and misses the spatial and temporal knowledge regarding the extent and dynamics of the surface water. The centroid layer of 266 North American Ramsar sites was used to identify the location of each site and determine the extent of each site based on the water occurrence results of this study. To demonstrate the capabilities of the implemented method and illustrate the details of water occurrence maps in 30 m spatial resolution, 33 sites were selected from wetlands and waterbodies. These sites were chosen to include different wetland classes, including inland, coastal, freshwater, saline, permanent, and temporary wetland types. The sites are located in different parts of North America, and 30 of them are designated as Ramsar sites. To show the results from the Caribbean, a Ramsar site from Cuba (Humedal Delta del Cauto) was included. Three sites (Ashepoo-Combahee-Edisto, Atchafalaya Delta, and Great Salt Lake) were selected from wetlands and water bodies that are not designated by the Ramsar Convention to depict the results of other regions as well. The maps in this study were produced for the entire North America since the Ramsar sites only represent the major and big wetland sites, and local smaller wetlands and other forms of surface water are not registered in this convention. For example, the few Ramsar sites in Canada do not cover a huge amount of wetland areas. Figure 12 shows the water occurrence map of 33 sites across North America.

Figure 12. Water occurrence map of selected wetlands and waterbodies in North America from 1982 to 2022. Details of the sites are provided in Table 2.

To demonstrate surface water dynamics, inundation class, and changes in the produced maps, results for Great Salt Lake are illustrated in Figure 13. The true-color summer median composite image of this region is shown, along with the DSWE classification results. The water occurrence map of this lake is produced in two timeframes, from 1982 to 2002 and from 2002 to 2022. Prior research highlights climate change and swift population expansion—Utah being among the fastest-growing and most arid states—as primary factors driving the situation (Baxter and Butler 2020). Throughout the twenty-year span of the western megadrought, water redirection from rivers supplying the lake has amplified to accommodate agricultural demands and the burgeoning urban areas’ increasing water consumption. This trend is illustrated in the results of the maps produced, as DSWE from 1985 to 2020 shows the shrinking of the lake, as well as water occurrence maps. Changes are highlighted in change detections, as well as showing the water loss area.

Figure 13. Surface water dynamic maps and change results of Great Salt Lake.

To quantify the results and changes in the surface water dynamic in wetlands and Ramsar sites, the total area of inundation class and change types are reported in Table 2. Due to the lack of an official borderline for Ramsar sites in North America, hydrological basins and the appropriate region in the water occurrence maps were used to determine the border of each site. Also, for coastal wetland sites, a buffer from the official coastline was considered. Permanently flooded areas show the area where the water was persistent most of the time, while temporarily flooded areas show the area where water happened only. Water gain and loss areas demonstrate the changes in wetland sites and depict the areas where water inundation frequency increased (wetter) or decreased (dryer). A wetland complex can experience both types of changes in different areas, however, in most wetland sites, only one change type is significant.

Table 2. Inundation class and changed area of wetland sites.

ID	Site Name	Country	Permanently Flooded (km^2)	Temporarily Flooded (km^2)	Water Gained (km^2)	Water Loss (km^2)
A1	Hay-Zama Lakes	Canada	177.76	567.95	4.64	85.27
A2	Whooping Crane Summer Range	Canada	13,657.31	1,566.26	115.70	213.35
A3	Peace-Athabasca Delta	Canada	10,413.07	1,896.19	27.70	757.62
A4	Dewey Soper Migratory Bird Sanctuary	Canada	15,296.60	1,002.13	388.92	586.87
A5	Fraser River Delta	Canada	2,128.32	426.80	6.98	6.57
B1	Creston Valley	Canada	94.78	8.39	4.56	3.71
B2	Beaverhill Lake	Canada	237.98	488.57	11.43	449.59
B3	Last Mountain Lake	Canada	215.65	48.67	46.85	0.94
B4	Quill Lakes	Canada	648.98	449.86	578.76	0.59
B5	Sand Lake National Wildlife Refuge	United States	88.21	80.24	101.61	11.84
C1	Humbug Marsh	United States	140.97	66.09	3.13	2.28
C2	Laguna de Santiaguillo	Mexico	120.33	64.60	36.21	4.31
C4	Lago de Chapala	Mexico	3,209.84	30.32	6.40	37.94
C5	Chesapeake Bay Estuarine Complex	United States	25,376.48	342.96	187.19	45.75
D1	Caddo Lake	United States	55.66	15.03	0.54	2.95
D2	Catahoula Lake	United States	96.70	81.40	9.87	20.47
D4	Atchafalaya Delta	United States	1,550.79	200.10	92.60	132.01
D5	Ashepoo-Combahee-Edisto	United States	55,999.49	851.01	260.25	539.79
E1	Everglades National Park	United States	22,347.89	2,825.98	672.93	1,511.90
E2	Humedal Delta del Cauto	Cuba	8,181.97	155.41	122.14	72.57
E3	Laguna San Ignacio	Mexico	3,833.58	83.86	58.76	29.24
E4	Laguna de Yuriria	Mexico	5,904.94	56.82	7.27	15.00
E5	Upper Mississippi River Floodplain Wetlands	United States	3,827.34	93.49	60.37	21.27
F1	Sistema Lagunar San Ignacio - Navachiste - Macapule	Mexico	15,224.93	251.73	126.01	85.70
F2	San Francisco Bay/Estuary (SFBE)	United States	77.26	43.52	69.10	0.62
F3	Marismas Nacionales	Mexico	11,766.99	232.61	163.63	56.95
F4	The Emiquon Complex	United States	1,171.74	85.27	208.91	16.70
F5	Laguna Madre	Mexico	201.01	146.30	31.79	47.48
G1	Ensenada de Pabellones	Mexico	32,501.64	902.80	534.96	135.73
G2	Reserva de la Biosfera Pantanos de Centla	Mexico	2,906.09	268.52	154.05	190.02
G3	Laguna Ojo de Liebre	Mexico	1,792.59	79.16	114.01	28.47
G4	Great Salt Lake	Mexico	7,205.95	1,885.57	180.71	2,786.32
G5	Humedales del Delta del Río Colorado	Mexico	11,508.04	601.86	583.81	39.58

Discussion

Leveraging the computational power of GEE and Landsat open-access historical archives for the past four decades, the 30 m surface water dynamic maps of North America were produced in this study. From coastal to inland wetland sites and different types of waterbodies, regardless of the temporary or permanent, saline or freshwater, spatial knowledge about extent, water occurrence, inundation class, and changes of wetland sites and waterbodies were provided in this study.

Landsat outlook

After 2008 and by removing the price of Landsat data and making it freely available for scientific and commercial usage, the true potential and value of this data were proven, and various applications and methods have been developed specifically based on Landsat data (Woodcock et al. 2008). Monitoring long-term changes requires not only the availability of historical observations but also the integrity and consistency of data quality through time and space. The proposed method in this study was based on calculation and processing for each pixel through a time series, though geometric, radiometric, and atmospheric corrections, along with robust cloud detection, are required to ensure a clear and consistent time series of observation for each pixel. Landsat Collection 2 Level-2 science-grade products provide the mentioned requirements with Tier 1 L1TP georeferencing standard with less than 12 m Root Mean Square Error (RMSE) that guarantee the consistent location of each pixel through the time series and minimize the errors and noises caused by mislocating the pixels. Considering the dependency of the DSWE unsupervised classification method on spectral information and to apply a model on a large geographic region with various atmospheric imaging conditions at different times, it is crucially important that a pixel correspondent to a certain type of land cover (i.e., water in this case) have a consistent value through time and space, and between different Landsat instruments and platforms. Landsat Level-2 surface reflectance data, as a physical value only dependent on the land cover type, regardless of atmospheric condition, Landsat instrument, and location, and is consistent within the time series. Radiometric and geometric consistency of Landsat data, along with a robust Fmask algorithm (Qiu et al. 2019) for cloud, cloud shadow, and snow detection, ensure a clear time series of observation on a pixel level to monitor dynamics of wetlands with the proposed method. Each Landsat platform collects images every 16 days, concluding between 22 and 23 images per year. However, due to Landsat policies, it is planned that at least two Landsat platforms should be operating, reducing the temporal coverage to 8 days. Currently, with the launch of Landsat-9 in 2021 and the continuity of Landsat-7 and -8, there is three available Landsat satellite that provides a suitable number of observations within a year for monitoring surface water dynamics and intra-annual changes in wetlands and waterbodies. However, due to different cloud coverage in different locations, producing annual hydroperiod maps of wetlands with Landsat data still remained a challenge. With the launch of Sentinel-2 by the European Space Agency and the availability of Harmonized Landsat Sentinel-2 (HLS) data (Claverie et al. 2018), the frequency of open-access science-grade optical data is greater than ever. This provides great opportunities for future studies with a focus on annual or short-term hydroperiod monitoring.

Methods and reference data

Various approaches for surface water monitoring were proposed in previous studies based on different types of earth observation data. For example, Synthetic Aperture Radar (SAR) data has the advantage of penetration through clouds and consistent frequency of data, in addition to the capability of monitoring water bodies beneath the canopies and vegetation as well as unique backscattering coefficient of water from other land covers (Amani et al. 2022). However, these advances are limited to monitoring the current wetlands and waterbodies status and short-term changes and are not capable of illustrating long-term changes and dynamic patterns of wetland surface water. Between the studies using Landsat data, the JRC Global Surface Water dataset provides complementary data regarding the transition of surface water (Pekel et al. 2016). However, compared to the proposed method in this study, the JRC dataset lacks confidence levels in water detection and provides only binary water/non-water maps (Soulard et al. 2020). In contrast, water areas are classified into four different classes in this study. Previous studies also did not have the advantage of Landsat-9 data, along with Landsat Collection-2 reprocessed data, with improved geometric and radiometric accuracy and maximum integrity of data between different Landsat platforms. Most previous studies only focused on surface water dynamics (Pickens et al. 2022, 2020); however, Ramsar wetland sites and other local sites, including inland wetlands, waterbodies, and coastal wetlands, were monitored in this study, and water occurrence maps were produced. In addition to the regional patterns of wetland surface water dynamics, site-based information extracted in this study can be used for local wetland conservation and decision-making in different countries in North America.

Collecting reference data has been a challenge in the remote sensing field, especially for large-scale studies. One approach is to collect ground truth data, which requires a lot of resources due to the inaccessibility of several regions and political barriers on a continental scale. In addition, collecting ground reference data demonstrates the current status of land cover only and lacks the historical aspect (Copass et al. 2018). In the case of surface water, historical data are usually available based on hydrological components of water bodies, such as depth, gage, and discharge data, and lack spatial knowledge. Collecting reference data based on high-resolution Google Earth images over North America is a time-consuming task that requires a team of experts for visual interpretation and is limited to the years 2000 and after. The absence and challenges of gathering reference data from non-remote sensing sources, especially in large-scale, long-term studies, add more limitations. To evaluate the results, previous studies are considered the only option to assess the agreement of the results.

Nature-based solution

The usage of nature as a solution to address various social and environmental challenges from a local to a global scale, such as urban heat islands, sea-level rise, and climate change, is usually referred to as an umbrella concept called Nature-based Solution (NbS) (Cohen-Shacham et al. 2016). Numerous ecosystem services provided by wetlands indicate their potential and importance to address a variety of social, economic, and environmental challenges. Biodiversity support, soil moisture and groundwater level regulation, coastal protection, carbon sequestration, flood regulation, and water quality protection are examples of natural services provided by wetlands (Thorslund et al. 2017). To highlight the best use of wetland ecosystems as a sustainable and cost-efficient nature-based solution for a range of current and future environmental and social challenges, spatial and temporal knowledge of these ecosystems is required. First and foremost, these valuable ecosystems should be conserved, degradation trends should be controlled, and restoration should begin. Surface water dynamics and change maps produced in this study provide a unique historical pattern of wetlands surface water and identify the location, magnitude, and class of changes, which are essential for detecting damaged wetland sites. Furthermore, for measuring the impacts of wetlands as a nature-based solution, it is required to analyze the correlation of wetlands’ spatial and temporal information with relevant environmental parameters, such as temperature and water quality. In this way, wetlands’ role as a nature-based solution can be proven, and the impact of wetland surface water loss and gain can be measured on the surrounding environment. For example, with the provided surface water dynamics information, the impact of coastal wetlands on shore protection and floods can be analyzed.

In the realm of remote sensing, the examination of interactions between wetlands and their surrounding environments, including vegetation, ice-covered areas, and permafrost, is of vital importance in understanding the intricate dynamics of these ecosystems. Wetlands often exhibit symbiotic relationships with nearby vegetation, as they provide essential hydrological and nutrient support, promoting biodiversity and ecosystem productivity (Hemati et al. 2021). Furthermore, wetlands situated in proximity to ice-covered areas and permafrost play a crucial role in regulating local and regional hydrological cycles, as they store and release water, influencing the thawing and freezing processes of adjacent frozen landscapes (Perreault et al. 2017). Remote sensing techniques enable us to observe and quantify these complex relationships over time and space, contributing to our understanding of wetland resilience and vulnerability in the face of environmental change and anthropogenic disturbances (Wulder et al. 2018). By employing advanced remote sensing methodologies, researchers can better assess the ecological implications of these interactions and support sustainable management and conservation strategies for wetlands and their surrounding ecosystems.

Future works

Surface water is an important indicator of wetlands, and by monitoring the dynamics and extent of surface water and waterbodies, wetland dynamics can be studied. However, wetland vegetation and surrounding land cover are considered important parts of wetland ecosystems and are responsible for most of the environmental services. Monitoring wetlands with the inclusion of multiple wetland classes besides water remained challenging due to the lack of suitable reference data for modeling and evaluating results through a time series. Bootstrapping (Champagne et al. 2014) and cross-validation (Friedl et al. 2010) are examples of alternatives for providing reference data. Deep learning methods are developing so fast, and recently advanced architectures, such as Transfer, Weakly-supervised, and Semi-supervised learning, require minimum or no amount of reference data (Hosseiny et al. 2024). However, the evaluation and uncertainty of models trained on data for the current situation and applied to historical data remain questionable. With the launch of Landsat-9 and the availability of Landsat-8 and Sentinel-2A and B, the temporal resolution of such data will be about 2.5 d (Wulder et al. 2019). In addition, SAR data are becoming more available; thus, there will be more opportunities for hydroperiod mapping using high-frequency optical data and data fusion approaches. Another promising area for future exploration would be studying the impact of precipitation, impervious surfaces, soils, and other possible explanatory variables on wetland changes. To depict the value and functionality of wetlands as Nature-based Solutions, water inundation frequency and hydroperiod maps can be direct empirical inputs of models, simulating flood events, urban heat islands, or sea-level rise.

Conclusion

Considering the importance and global value of wetlands, there is limited spatial knowledge about wetlands surface water dynamics and changes. The surface water dynamic maps of North America were produced in this study, including water occurrence, inundation class, change magnitude, and change class. Leveraging the historical archive of the Landsat program and the cloud computational power of Google Earth Engine, 40 years of surface water dynamics in waterbodies and 266 Ramsar wetland sites were monitored with 30 m spatial resolution. Science-grade products of Landsat with superior radiometric and geometric consistency and cloud detection algorithms ensure the quality of time series to minimize the noises in time and space dimensions. The capabilities of unsupervised classification of the Dynamic Surface Water Extent method for water detection in the North American region were illustrated in this study, and the dynamics of wetlands were monitored using this method between 1982 and 2022. Using this method, water classes were detected in different confidence levels and achieved an average overall agreement of 92.97% with the Pekel et al. (2016) results. In addition to the inundation class and extent of wetland surface water and waterbodies, the dynamics and associated changes were detected by calculating water occurrence in the time series. The collection of surface water products and dynamic maps produced in this study are beneficial for wetland conservation programs and governments. Using the dynamics of wetlands surface water and the change detection provided in this study, the impact of wetlands as a nature-based solution for environmental changes can be measured. However, the use of multi-sensors for producing hydroperiod maps and the lack of ground-based reference data for applying machine/deep learning algorithms are still considered the main challenges in this field.

Acknowledgments

The authors would like to thank the United States Geological Survey (USGS) and Google for providing access to the Landsat data and the Earth Engine platform, which were instrumental to the success of this research.

Disclosure statement

The authors declare no conflicts of interest in relation to the research.

Data availability statement

The open-access datasets used in this study are freely available on various platforms and websites. The sources of the primary datasets used in the study are as follows:

• Landsat Collection 2 Surface Reflectance dataset: https://developers.google.com/earth-engine/datasets/catalog/landsat (Accessed 28 November 2023).

• ALOS World 3D – 30m (AW3D30) digital surface model (DSM) dataset: https://developers.google.com/earth-engine/datasets/catalog/JAXA_ALOS_AW3D30_V3_2 (Accessed 28 November 2023).

• Ramsar Sites Information Service wetland dataset: https://rsis.ramsar.org/ (Accessed 28 November 2023).

• Products produced in this study can be accessed through this link: https://github.com/BEEILAB/Publications/tree/main/Comprehensive%20Landsat-based%20Analysis%20of%20Long-term%20Surface%20Water%20Dynamics%20over%20Wetlands%20and%20Waterbodies%20in%20North%20America.

Additional information

Funding

The Authors acknowledge the Memorial University of Newfoundland and the research funding provided by VPR/SGS Pilot program, the Natural Sciences and Engineering Research Council (NSERC) Discovery (Grant No. RGPIN-2022-04766).