Forty Years of Wetland Status and Trends Analyses in the Great Lakes Using Landsat Archive Imagery and Google Earth Engine

: Wetlands provide many beneﬁts, such as water storage, ﬂood control, transformation and retention of chemicals, and habitat for many species of plants and animals. The ongoing degradation of wetlands in the Great Lakes basin has been caused by a number of factors, including climate change, urbanization, and agriculture. Mapping and monitoring wetlands across such large spatial and temporal scales have proved challenging; however, recent advancements in the accessibility and processing efﬁciency of remotely sensed imagery have facilitated these applications. In this study, the historical Landsat archive was ﬁrst employed in Google Earth Engine (GEE) to classify wetlands (i.e., Bog, Fen, Swamp, Marsh) and non-wetlands (i.e., Open Water, Barren, Forest, Grassland/Shrubland, Cropland) throughout the entire Great Lakes basin over the past four decades. To this end, an object-based supervised Random Forest (RF) model was developed. All of the produced wetland maps had overall accuracies exceeding 84%, indicating the high capability of the developed classiﬁcation model for wetland mapping. Changes in wetlands were subsequently assessed for 17 time intervals. It was observed that approximately 16% of the study area has changed since 1984, with the highest increase occurring in the Cropland class and the highest decrease occurring in the Forest and Marsh classes. Forest mostly transitioned to Fen, but was also observed to transition to Cropland, Marsh, and Swamp. A considerable amount of the Marsh class was also converted into Cropland.


Introduction
Wetlands provide many services to the environment and humanity, including water purification, protection from natural hazards, soil and water conservation, and shoreline protection [1]. Moreover, wetlands are capable of significant water storage, serve as a form of natural flood control [1], transform and retain chemicals [2,3], and provide habitat for many species of plants and animals [4].
Wetlands are degrading in the world at a fast pace [5,6]. These valuable natural resources are being affected by urbanization, irrigation, climate, pollution, and other anthropogenic activities [5,7,8]. Thus, their ecological functions are negatively influenced, resulting in habitat loss and biodiversity reduction [8]. For instance, Ref. [9] reported that more than 50% of wetlands in North America and Europe are degraded. In Canada, this trend is similarly observed and is expected to continue in the future [10]. However, it is challenging to accurately estimate the extent and type of wetland loss using traditional ground-based approaches. An alternative approach to assessing large-scale wetland changes over time is through remote sensing methods [5,6].
Numerous remote sensing systems capture consistent repetitive imagery with global coverage. Satellite datasets have frequently been used to map wetlands across Canada and have proved to have a high potential to delineate various wetland types [7,[11][12][13][14][15][16][17], as well as monitor wetland changes [5,6,18]. For example, optical satellite imagery is commonly used for wetland mapping due to its ability to distinguish spectral differences in wetland species [19]. In this regard, Landsat imagery in particular also provides a historical record of earth observation data from 1972 to the present, allowing for the monitoring of wetland changes occurring over past decades, particularly where field data have not been consistently collected.
Satellite imagery has become increasingly accessible, along with improved processing efficiency. Until recently, acquiring imagery across large spatial and temporal scales was slow and challenging. However, recent advancements have efficiently allowed the acquisition of multiple sources of imagery for large-scale wetland mapping and monitoring initiatives. For instance, Google Earth Engine (GEE) is a cloud-based platform which provides access to a wide range of earth observation products, including the Landsat historical record. In addition to the availability of pre-processed imagery, several machine learning techniques have been implemented and can be imported and modified by users [20][21][22]. To date, many researchers have used GEE to effectively analyze wetland changes over long periods of time. For example, Ref. [6] mapped wetland changes in Newfoundland, Canada, over a 30-year period using the Landsat historical archive in GEE [6]. Additionally, Ref. [5] employed the Landsat archive to develop a machine learning model in GEE that classified wetlands with approximately >80% accuracy for each class [5]. They also identified the percentage of wetland areas that had undergone changes since 1984 for the province of Alberta, Canada.
The Great Lakes region and its subbasins encompass approximately 765,000 km 2 across the Canada and United States border and contains important wetland areas. Although the region provides millions of people with fresh drinking water, food, and recreational opportunities, it is highly impacted by human disturbances and subsequent wetland loss. Although multiple wetland mapping initiatives have already been undertaken in the Great Lakes basin (e.g., [23,24]), there has not been a comprehensive assessment of wetland change over the entire basin. For example, a user-friendly interface, named Wetland Extent Tool (WET) was developed within GEE by [23] to classify wetlands in the Great Lakes basin. WET was a semiautomated wetland classification workflow to classify multisource data (i.e., optical satellite, SAR, and elevation images) using the Random Forest (RF) algorithm. WET used a pixel-based classification approach and, thus, was affected by noise due to within-class spectral variation. Moreover, it only resulted in three general categories of wetland, upland, and open-water, and did not include wetland subclasses (i.e., Bog, Fen, Marsh, and Swamp) which are common in Canada and the US. In another study [24], dynamic changes in wetlands, surface water extent, and water level change of wetlands in a portion of the Great Lakes were monitored between 2016 to 2018 using SAR, multispectral, elevation, and InSAR products. The authors mapped dynamic surface water and flooded vegetation extent by thresholding the Radarsat-2 HV polarization observations. The water level monitoring was also conducted using InSAR analysis with the C-band HH polarization SAR data. Additionally, they employed a variety of approaches for the classification and change detection of wetlands in the Great Lakes basin, including spatiotemporal classification, Object-Based Image Analysis (OBIA), and SAR-based classification. Although the authors of [24] exhibited promising results for wetland change detection in the Great Lakes basin, it was conducted over a few selected pilot sites and did not extend to the entire basin or multiple decades. Furthermore, such an approach requires high-resolution datasets, which are usually expensive and not freely available.
In this study, we developed a model to examine changes in wetland types and extents in the Great Lakes basin over the past four decades. To this end, advanced classification and Change Detection (CD) algorithms were developed and applied to the archived Landsat imagery within the GEE cloud computing platform. In contrast to the authors of [23], an object-based RF classification was used to generate a more accurate wetland map. Moreover, four wetland classes (i.e., Bog, Fen, Marsh, and Swamp) and five non-wetland classes (i.e., Open Water, Forest, Grassland/Shrubland, Cropland, and Barren) were delineated to provide more detailed wetland maps. Finally, open-access Landsat data were employed to facilitate wetland change analysis over large areas and multi-decade time frames.

Study Area
The Great Lakes ( Figure 1) are a series of large and interconnected freshwater lakes in central North America, whose individual subbasins form the larger Great Lakes basin and constitute the largest freshwater reserve on Earth. The Great Lakes basin was formed during the recession of the Laurentide ice sheet and now contains diverse ecosystems and abundant plant and wildlife [25,26]. The region spans a large and variable area, but generally hosts a humid continental climate, which can be influenced by air masses from other regions (e.g., cold Arctic air, warm tropical air). Wetlands occur throughout the study area, most commonly in coastal and riverine settings, though they can also be found inland.  Figure 1. The samples were well distributed over the Great Lakes basin and, thus, these samples can properly represent different types of wetland and non-wetland classes in the study areas. Table 1 provides information regarding the number and area of the field samples collected and prepared as part of this study. In total, nine landcover classes were considered based on the Canadian Wetland Classification System (CWCS) system and previous wetland studies conducted in Canada and the US: four wetland classes of Bog, Fen, Marsh, and Swamp, as well as five non-wetland classes of Open Water, Forest (e.g., deciduous, coniferous, and mixed woodlands), Grassland/Shrubland, Cropland, and Barren (e.g., urban, rock, bare soil, sand, and other non-vegetated areas). The majority of sample areas belonged to the Marsh and Forest classes, respectively, while the number and area of samples from the Bog and Grassland/Shrubland classes were relatively small. It is worth noting that the Open Water class was composed of both Shallow Water and Deep Water classes. However, due to the low number of Shallow Water field samples and their lower accuracy, it was decided to consider one water class as Open Water. The study area (the Great Lakes basin and five subbasins of Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan) and the distribution of field data. The red box shows the location of the study area. The distribution of all field samples is illustrated in Figure 1. The samples were well distributed over the Great Lakes basin and, thus, these samples can properly represent different types of wetland and non-wetland classes in the study areas. Table 1 provides information regarding the number and area of the field samples collected and prepared as part of this study. In total, nine landcover classes were considered based on the Canadian Wetland Classification System (CWCS) system and previous wetland studies conducted in Canada and the US: four wetland classes of Bog, Fen, Marsh, and Swamp, as well as five non-wetland classes of Open Water, Forest (e.g., Figure 1. The study area (the Great Lakes basin and five subbasins of Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan) and the distribution of field data. The red box shows the location of the study area.

Satellite Data
In this study, the archived level 2 products (surface reflectance) of Landsat-5, -7, and -8 images from 1984 to 2021 available in GEE (https://developers.google.com/earthengine/datasets/catalog/landsat, accessed on 20 January 2022) were employed ( Table 2). In total, 58,509 Landsat images were acquired and processed. As outlined in Table 2, the images from every three-year period (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998) or every two-year period  were combined to produce cloud-free Landsat mosaic images for the entire study area. It should be noted that the study area is covered by clouds and snow for much of the year and creating cloud-free images with less time intervals (T), such as annual images, was not possible. In total, 17 time intervals were considered and, thus, 17 wetland maps were produced to analyze changes over the past four decades.

Methodology
The employed methodology included four main steps. First, the field data were preprocessed to make them suitable for wetland mapping and change analysis (see Section 2.4.1). Subsequently, all available field samples were analyzed in order to select only spectrally unchanged samples for use in the classifications of all time intervals (see Section 2.4.2). Then, the selected samples were used to train the classification algorithm and produce wetland maps at different time intervals (see Section 2.4.3). Finally, a CD algorithm was applied to assess wetland trends over the past four decades (see Section 2.4.4).

Field Data Preparation
The point-based samples were used to delineate the corresponding site boundaries using high-resolution satellite images available in Google Earth and ArcGIS (see Figure 2 for an example). Delineated polygons were subsequently reassessed for accuracy (in terms of both location and extent) before being included in the classification model.

Unchanged Field Samples Selection
The flowchart of the method to generate spectrally unchanged field samples is illustrated in Figure 3. In this study, the Continuous Change Detection and Classification (CCDC; [27]) method was utilized to select unchanged samples and generate additional field samples for each time interval when field data had not been collected during those years. The CCDC method analyzes all field samples and selects the samples whose spectral responses have remained unchanged over the study period; in this case, four decades. The resultant samples are referred to as the unchanged field samples and used as training and test data for all years. First, cloud, cloud shadow, and snow/ice pixels were masked in the Landsat images using the quality band of Landsat (i.e., Bitmask for pixel_qa) and C Function of Mask (CFMASK) algorithm [28] within GEE. The quality band flags undesired pixels (e.g., cloud and snow) with a class name or unused designation, depending on the Landsat generation [29]. The masked Landsat images were then applied to calculate three spectral indices, Normalized Difference Vegetation Index (NDVI, Equation (1)), Normalized Difference Water Index (NDWI, Equation (2)), and Normalized Difference Build-up Index (NDBI, Equation (3)), which were then used to determine the unchanged field samples using the CCDC model derived in Equation (4).
where i, x, and t, respectively, indicate the spectral index, Julian date, and the number of days per year (i.e., 365.25 days); a 0,i stands for the overall value of the spectral index i of a Landsat image; a 1,i and b 1,i specify the intra-annual change; and c 1,i shows the interannual change. The model coefficients were determined from the field samples. The new samples were inserted into Equation (4) to estimate model values and residuals by comparing the observed and modelled sample values. A threshold value was set to 20%, where it was assumed that an abrupt or interannual change had occurred if residuals exceeded 20%. Samples that did not meet the threshold were not used in the classification of another time interval. The final remaining samples were assumed to have stable spectral responses and therefore had not undergone change over the study period. As such, these samples were used to classify wetlands for all 17 intervals. Prior to classification, the unchanged field samples were randomly split into training (70%) and test (30%) groups, where the training group was used to train the classification algorithm and the test group was used to assess the accuracy of the resultant wetland maps.

Classification
The steps of the classification method to produce wetland maps at different time intervals are demonstrated in Figure 4.
Wetlands maps for each time interval were generated using the unchanged field samples identified through the CCDC process. Cloud, cloud shadow, and snow/ice were masked from all images, which were then divided into two groups based on season: (1) the months of April, May, June, and July, which was called the Spring-Summer time, and (2) the months of August, September, and October, which was called the Summer-Fall time. The reason for this was because it has been widely reported that multi-temporal datasets could improve wetland detection [7,13,30]. All images within each group were then down-sampled into a single mosaic by calculating the mean value per pixel over the entire time series (i.e., the Spring-Summer and Summer-Fall mosaic images). Ten features (layers), which included the seven main spectral bands of Landsat images (i.e., blue, green, red, Near Infrared (NIR), two Shortwave Infrared (SWIR), and Thermal Infrared (TIR) bands), and the three spectral indices (NDVI, NDWI, and NDBI) were extracted from each mosaic image and used in the classification. Object-based classification techniques have been widely reported to produce better results when compared to pixel-based methods [13,30,31]. Additionally, large heterogeneous landcover types, such as wetlands, have been shown to require an object-based approach [11]. Therefore, the mosaicked images were ingested into the Simple Non-Iterative Clustering (SNIC) algorithm to be segmented. The SNIC method evenly distributes a number of seeds throughout the image to partition it into super-pixels. A priority queue is then applied to determine the next pixel to be assigned to a cluster based on the distance of the pixel from the segment centroid. After each new pixel is added to a cluster, a new centroid is computed. This process continues until centroid convergence [32].
The final segmented image, consisting of 20 layers (i.e., 10 layers from each of the Spring-Summer and Summer-Fall mosaic images), was ingested into an RF classification algorithm. RF has proven to be an effective classifier for landcover classification, especially for wetland mapping [7,13,30,33]. RF contains a set of decision trees that divide the input pixels into mutually exclusive groups until each node represents one of the final classes. Model tuning parameters are selected based on the available samples and the objectives of the classification. In this study, these parameters were selected after several trial-and-error classification iterations and based on the results of previous wetland mapping studies within GEE [5,34,35]. Accordingly, the following values were selected for each tuning parameters of the RF algorithm: numberOfTrees = 50, variablesPerSplit = 5, minLeafPopulation = 1, and bagFraction = 0.5.

Accuracy assessment
After the classification was conducted, it was visually assessed using Google Earth high-resolution images and basemaps available in ArcGIS. If the visual inspection suggested the classification could be further improved, the classification parameters and/or the training samples were revised, and the RF retrained until satisfactory results were obtained. For the statistical accuracy assessment, the final classification was compared with the test data and the results presented in a confusion matrix. The confusion matrix was used to examine whether different classes in the generated wetland maps were confused with other classes.

Change Detection (CD)
Change analysis was conducted using the 17 produced wetland maps, where the map produced for 2020-2021 (T17 in Table 2) was selected to be the reference time interval as it covered the most recent time interval and achieved the highest map accuracies. It should be noted that the selection of the reference interval does not significantly affect the analysis and is simply a starting point for the CD algorithm.

• Change Detection (CD) between two time intervals
The differences between every two time intervals were first determined using two CD methods of class differencing and image differencing (see Figure 5). The combination of the two methods improved CD results as opposed to only using a single method. Furthermore, utilizing pixel-and object-based CD methods decreased the overestimation of changed areas and noise, respectively.
The pixel-based change analysis was conducted using the spectral bands of the Spring-Summer and Summer-Fall mosaic images, where the spectral distance between the mosaicked images of two time intervals was calculated. A threshold of 70% was subsequently applied to the resultant map, where a difference greater than 70% between the spectral responses of pixels at two different time intervals indicated changed pixels. The object-based change analysis was conducted by differencing two object-based wetland maps. The final change map was produced as the intersection of the pixel-and object-based change maps, where the final change map was object-based.

• Change Detection (CD) through all time intervals
As the main objective of this study was to analyze the changes occurring over the past four decades, an additional analysis was conducted (see Figure 6) to determine overall wetland changes occurring during the past 40 years. First, using the reference map (i.e., T17) as one of the two change analysis maps, the other map (i.e., T16) was updated. Specifically, the class labels of the changed objects in the other map were updated based on the reference map class labels. Through this approach (Figure 6), all 16 other maps were hierarchically updated, and the overall changes were determined.  Table 2).  Table 2) and CD is the method proposed in Figure 5.  Figure 7 presents the generated wetland map of the reference time interval (i.e., 2020-2021) for the entire Great Lakes basin and its five subbasins individually.

Results and Discussion
Visual inspection of the produced wetland maps showed satisfactory accuracy in identifying different wetland and non-wetland classes. For example, the method accurately delineated large water bodies and urban areas as Open Water and Barren. Additionally, the southern portions of the Great Lakes basin, which are dominated by agricultural lands, were correctly identified as Cropland. The quantitative values extracted from the reference wetland map indicated that the study area is mainly covered by Open Water (34.8%), Forest (24.6%), and Cropland (17.7%), respectively. On the other hand, Bog (0.003%), Grassland/Shrubland (0.7%), and Swamp (1.25%) were the classes that covered the least amount of area within the Great Lakes basin. Further visual investigations at subbasin levels also suggested that the methodology was capable of capturing smaller water bodies and urban areas.
Visual intercomparisons using high-resolution satellite images were also carried out to assess the quality of the produced wetland maps at smaller scales. Two zoomed-in classified areas, along with their corresponding high-resolution satellite imagery, are demonstrated in Figure 8. Based on Figure 8a, the implemented approach correctly distinguished different land covers. For example, the agricultural lands outside and surrounding metropolitan areas were correctly identified as Cropland. Likewise, the urban areas, water bodies, and forest patches were also accurately assigned to Barren, Open Water, and Forest, respectively (see Figure 8b). Considering both zoomed-in patches, several pixels located around large water bodies and multiple small shallow water areas were also correctly classified as Marsh (i.e., emergent marshes), which was also observed in other parts of the produced wetland maps.
Despite the high visual accuracy of the maps, there were errors, such as the misclassification of Forest as Swamp and incorrect classifications of Cropland as Grassland/Shrubland or other wetland classes. The confusion between the Forest and Swamp classes has been reported in many other wetland mapping studies [5,36]. Swamp and Forest have similar spectral responses in optical satellite imagery; thus, it is challenging to separate these two classes when only optical data are employed [36]. Other remote sensing datasets, such as L-band Synthetic Aperture Radar (SAR) data, are required to better discriminate between these two classes. Figure 7. Generated wetland map of the (a) Great Lakes and (b-f) its subbasins (i.e., Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan, respectively) using satellite imagery between 2020 and 2021 as the reference time interval.

Accuracy Levels
In addition to the visual evaluation of all wetland maps, independent test samples (30% of field data) were employed to quantify the accuracy of each generated wetland map by calculating a confusion matrix. As illustrated in Figure 9, all of the wetland maps had OAs and KCs greater than 84% and 81%, respectively, indicating the high potential of the proposed methodology for large-scale wetland mapping using only spectral indices derived from optical satellite data. The low variation in the obtained OAs (85.2% ± 0.7%) and KCs (82.4% ± 0.8%) also demonstrated the consistent performance and robustness of the applied classifier. As is clear in Figure 10, the average UAs of all wetland maps were between 76.6% and 80.3%, with the lowest and highest average UAs belonging to wetland maps of the 13th (i.e., 2013-2014) and 8th (i.e., 2003-2004) time intervals, respectively. In terms of the class-wise UAs, Open Water and Cropland obtained the two highest UAs in all time intervals with average UAs of 99.1% and 95.8%, respectively, followed by Forest, Barren, Marsh, Fen, Swamp, Grassland/Shrubland, and Bog. Bog acquired the lowest UAs in almost all time intervals, except for the 4th time interval, with UAs between 48% and 57%. However, Bog had higher PAs in most time intervals in comparison to other classes, leading to average PAs of approximately 99.6%. Open Water and Grassland/Shrubland had the second-and third-highest PAs with average values of 99.2% and 97.3%, respectively. The lowest PAs were recorded for Marsh in almost most of the generated wetland maps, with an average PA of 79.1%. Overall, the average PAs were relatively higher than the average UAs, commencing at 87.6% in the 13th (i.e., 2013-2014) time interval and peaking at 90.1% in the 1st (i.e., 1984-1986) time interval. The confusion matrices for all wetland maps were also calculated. As an example, Table 3 provides the confusion matrix of the wetland map for the last time interval (i.e., 2020-2021). Overall, the largest confusion was observed between the wetland classes (i.e., Bog, Fen, Marsh, Swamp), which was anticipated due to their comparatively higher spectral similarities. This high spectral similarity is primarily associated with their similar ecological conditions [19]. Among wetland classes, the highest confusions were observed for the Fen class, in which 108 pixels (9.8% of all field samples of Fen) were incorrectly classified as other wetlands. Marsh had the second-most considerable confusion with other wetland classes (~8%). The confusion matrix also illustrated that the Open Water class had the lowest misclassification rates among other classes, owing to its distinct spectral characteristics. There were also confusions between wetland and non-wetland classes, the highest of which was observed between Marsh and Barren, Cropland, and Grassland/Shrubland. This issue was also noticed in reverse, causing the Marsh class to have a low PA and UA simultaneously in comparison to other classes. Moreover, Bog was confused with other classes, in which multiple pixels (77 out of 151 pixels) of other classes were misclassified as Bog, leading to the lowest UA value for this class. This pattern was also observed for the Grassland/Shrubland class with the second-lowest UA value, in which almost 42.1% of pixels from other classes were incorrectly classified as Grassland/Shrubland. Furthermore, there was confusion across non-wetland classes, with the most (130 pixels out of 3845) occurring between Barren and Cropland.

Change Analysis
All of the generated wetland maps were intercompared and several CD analyses were performed to explore wetland changes between 1984 and 2021. Since the generated wetland maps had acceptable accuracies, they were assumed to provide reliable results during the CD analyses. However, the intrinsic misclassification errors, such as confusions between different classes that reduced the PAs and UAs of the classes, could introduce uncertainties in the CD results (see Figure 10). For instance, Marsh, Forest, Barren, and Fen had relatively lower PAs than other classes (Figure 10b), which could be due to a high similarity between the spectral signatures of these classes. Moreover, the temporal changes of wetlands, such as seasonal changes, had considerable impacts on the results of both classification and change analysis.

Overall Change
As part of the initial CD analysis, all wetland maps were intercompared through the post-classification procedure to produce a binary map depicting changed and unchanged areas within the Great Lakes basin (see Figure 11). Visually, the Open Water portions of the Great Lakes basin were not subjected to change, while a regular change pattern was not seen in the study area. However, several hotspot areas of alteration were observed over the north/northwest, west, east/south, east, and west regions of the Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan subbasins, respectively. These findings revealed that changed areas constituted approximately 16.39% (125,101 km 2 ) of the Great Lakes basin, while 83.61% of the study area (638,125 km 2 ) remained unchanged. Further investigations at the subbasin level suggested that the highest changed area was in the Lake Ontario subbasin, whereas the least amount of change was in the Lake Erie subbasin, with respect to the total area of each subbasin. In particular, the Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan subbasins were subjected to 15.15%, 14.71%, 23.65%, 14.14%, and 17.88% changes over the past four decades, respectively. Visual inspections also revealed that several changed areas were related to urban and agricultural expansion, as well as deforestation, probably due to anthropogenic and human-induced activities. However, as discussed, the errors of the classified maps resulted in some uncertainties in the change results. For example, multiple false detection and noisy areas in the final change maps, which especially happened between forest and wetland types, were the results of errors in the produced wetland maps. Another reason for potential errors in the change maps could be the effects of the seasonal changes in wetland classes.

Change Trend Analysis
The computed class areas from the generated wetland maps were examined to quantify the relative changes in the area of each class. To this end, the proportionate area of each class in the earliest time interval was considered as the starting point, and the relative differences were calculated (see Figure 12). The results demonstrated that Open Water and Bog, the majority and minority land covers of the Great Lakes basin, respectively, covered the most consistent amount of area during the study period. Similarly, Barren, Grassland/Shrubland, Fen, and Swamp had relatively stable areas between 1984 and 2021, with slight upward trends for the first three aforementioned classes and a slight downward trend for the last class. The highest increase corresponded to Cropland, with an approximately 3.2% growth in area, indicating the significant impact of anthropogenic practices in this region. In contrast, Marsh and Forest had the two highest downward trends with almost 2.8% and 2.1% losses, respectively, over the past four decades. Overall, the increase and decrease in the areas of Cropland and Forest were expected due to several anthropogenic activities, including agricultural expansion and deforestation, that occurred in the Great Lakes region, especially over recent years. Figure 11. The produced changed and unchanged binary map of the (a) Great Lakes and (b-f) its subbasins (i.e., Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan, respectively) between 1984 and 2021.
The areas of each class, derived from the 17 wetland maps, are also presented in Figure 13 to better quantify each class's dynamic between 1984 and 2021. Figure 13 shows that Open Water, Forest, and Cropland were the most dominant classes in all time intervals, covering areas of approximately 265,500 km 2 , 195,300 km 2 , and 128,300 km 2 , respectively. Regarding wetland classes, the Great Lakes basin was covered mainly by Marsh and Fen, with areas of approximately 78,100 km 2 and 56,200 km 2 , respectively. Overall, Barren, Cropland, Fen, and Grassland/Shrubland presented steady rising trends with varying slopes, while the Forest, Marsh, and Swamp classes had stable downward trends. For instance, Barren areas started at almost 23,400 km 2 in the first time interval and increased to more than 26,400 km 2 in the latest time interval, which is mainly linked to urban growth in metropolitan areas such as Toronto, Buffalo, and Green Bay. Additionally, the Cropland areas exceeded 135,200 km 2 , with an almost 12.2% increase during the study period, and could be considered another human-induced indicator for wetland changes in the study area. Fen experienced a massive increase among wetland classes, almost 10,000 km 2 , reaching an area greater than 61,600 km 2 . Conversely, Marsh experienced the highest rate of decrease, falling from 92,600 km 2 to 70,700 km 2 . Likewise, the Forest area declined from almost 203,900 km 2 to approximately 188,100 km 2 in the past four decades, with anthropogenic activities such as urban and agricultural expansion, as well as misclassification errors in the wetland maps, the likely cause. Moreover, the significant number of forest wildfires recorded in the past few decades could be another contributing factor in the loss of forest area. In particular, the Canadian National Forestry Database and the United States Environmental Protection Agency reported many wildfires events in Ontario, Canada, and Michigan, US, that resulted in notable forest losses [37]. Meanwhile, Open Water and Bog, two classes with minor relative changes (see Figure 12), had fluctuating trends with various downward and upward trajectories, culminating in an increase and a decrease in their area at the end of the study period, respectively.

Gain and Loss
The generated binary map of the changed and unchanged areas ( Figure 11) was manipulated to illustrate the spatial pattern of wetland gains and losses ( Figure 14). The Gain (i.e., wetland gain) is associated with non-wetland classes that were converted to either one of the wetland classes, whereas the Loss (i.e., wetland loss) indicated wetlands that were converted into non-wetland classes. In total, the Great Lakes basin was subjected to 6.8% wetland loss and 4.9% wetland gain between 1984 and 2021. According to the Gain and Loss maps (Figure 14), the Lake Superior subbasin had a higher wetland gain (6.2%) than wetland loss (4.8%) during the past four decades. However, wetland losses exceeded wetland gains in the four other subbasins, with approximate differences of 0.6%, 6.0%, 4.8, and 3.8% for Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan, respectively. Moreover, the loss of coastal Marsh and Swamp were observed in several parts of the Great Lakes basin between 2016 to 2018 ( Figure 14), which was in line with the findings of the precious research study [24]. Additionally, coastal wetland loss in some parts of the Great Lakes basin (e.g., Lake Erie and Lake Michigan) were in good agreement with the findings of previous studies [38].
In a detailed quantitative overview, the gain and loss values were also computed for all classes and are illustrated in Figure 15. The Marsh and Forest classes had the two greatest losses of 44,408 km 2 and 42,406 km 2 , while they also experienced gains of 23,227 km 2 and 26,866 km 2 , leading a net loss of 21,181 km 2 and 15,540 km 2 over the past four decades, respectively. Swamp area also experienced a net decline, losing 8495 km 2 and gaining only 4950 km 2 , making it the second-largest source of wetland loss in the Great Lakes basin. Meanwhile, Bog, Open Water, and Grassland/Shrubland had the lowest gain and loss values, respectively, which explains the observed relative stability of these classes over the past four decades. Cropland experienced major gains (29,896 km 2 ) and a minor loss (5279 km 2 ), resulting in a net increase of 24,617 km 2 . Figures 16 and 17 illustrate the wetland gain and loss maps at different time intervals, respectively. According to Figure 16a   ) its subbasins (i.e., Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan, respectively) in the past four decades. According to Figure 17a, the most significant loss in the Great Lakes basin occurred between 1996 and 1998, closely followed by the loss rate between 2013 and 2014. However, the least amount of loss occurred between 2020 and 2021. In Lake Superior, the greatest loss occurred during the 1990-1992 interval, while the smallest loss occurred during the 2019-2021 interval (Figure 17b). It was also observed that most of the Lake Huron

Transition between Classes
A From/To map (Table 4) was also generated to explore the transition rate between each pair of classes. To this end, the first and the last wetland maps of the Great Lakes basin were compared to determine the amount of land cover transition that occurred between the nine classes. The most significant transition was observed between Forest and Fen, where approximately 18,449 km 2 of Forest in 1984-1986 was converted to Fen by 2020-2021. Both Open Water and Bog did not experience significant transitions from or to other classes. However, several classes experienced multiple land cover transitions over the past four decades. For example, approximately 17,765 km 2 , 14,158 km 2 , and 6965 km 2 of Marsh were converted to Cropland, Forest, and Fen during the past four decades, respectively. Among other wetland classes, Fen and Swamp were both primarily transformed into Forest or Marsh. Overall, wetland classes were most commonly converted to Forest or Cropland. The two main land cover origins for an increase in Barren area were found to be Marsh and Cropland, where 2496 km 2 and 4624 km 2 of each were converted to Barren, respectively. Furthermore, it was observed that an area of 22,480 km 2 was converted to Cropland between 1984 and 2021, mainly from Marsh and Forest, highlighting the agricultural expansion occurring in this region. Figure 16. Wetland gain maps of the (a) Great Lakes basin and (b-f) its subbasins (i.e., Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan, respectively) at different time intervals. Figure 17. Wetland loss maps of the (a) Great Lakes basin and (b-f) its subbasins (i.e., Lake Superior, Lake Huron, Lake Ontario, Lake Erie, and Lake Michigan, respectively) at different time intervals.

Limitations and Suggestions
In this study, due to the low number and quality of the samples, the Shallow Water class was combined with the Deep Water class to form the Open Water class. However, to conform the five wetland classes specified by the Canadian Wetland Classification System (CWCS, [39]), Bog, Fen, Marsh, Swamp, and Shallow Water, it is suggested additional Shallow Water field samples be collected in order to properly discriminate this class from Deep Water in future studies.
As explained in Section 2.4.2, the CCDC method was used to select the spectrally unchanged samples over the study period. During this process, it was observed that there were a considerable number of samples that had spectrally changed over the past four decades and, thus, were not used in the training and validation of the classification model. This reduced the number of samples for some classes, notably Bog, Grassland/Shrubland, and Swamp, and, consequently, affected the mapping accuracy of these classes. Therefore, similar to the Shallow Water class, it is suggested that more field data be collected and processed for these classes to improve the accuracy of the maps and change detection results. An alternative solution would be using wetland field data collected from other provinces/states of Canada/US with similar landscapes to the Great Lakes region.
The Landsat archive (i.e., Landsat 5, 7, and 8) was employed in this study to map wetlands and assess changes, as it is the only remote sensing system which provides continuous earth observation imagery over the four-decade period. However, Landsat imagery has a spatial resolution of 30 m and, as a result, introduces a source of misclassification in the change analysis. Moreover, discriminating multiple wetland classes using only optical imagery is challenging. For example, SAR data provides additional information to improve the separation of similar classes, such as Swamp and Forest [16,18,40]. Unfortunately, historical SAR data spanning the study period is not available and could not be applied to all classifications. Furthermore, since multiple images were required to produce a cloud-free mosaic image for each time interval (see Table 2), variations between different images could introduce another source of error.
Forest changes may result from clear-cutting or wildfires that mistakenly appeared as changes. It is recommended to further investigate the burned areas and/or clear-cut areas to obtain a more accurate understanding of reasons for change, especially in vegetated areas. To this end, vegetation change detection methods (e.g., [41]) along with our proposed approach would be beneficial.
Wetlands are naturally complex ecosystems that vary between and within individual sites due to changes in vegetation, water, and chemicals related to seasons, weather, and disturbance. The dynamic interaction between wetlands and their surrounding environments could also result in wetland features appearing different between years, although no change has in fact occurred [1].
In summary, there are multiple solutions to improve the accuracy of the results reported in this study. For example, the integration of high-resolution satellite images and accurate elevation information could considerably improve the accuracy of the classification and, subsequently, change the analysis. On the other hand, adding textural and geometrical features to the classification process could help to better separate different wetland types. Moreover, the results of previous studies [23,24,42] showed that integration of multi-temporal SAR and optical images in classification and change detection of wetlands resulted in greater accuracy. The Great Lakes undergo distinct water level variations that affect the water extend [38,43]. Thus, multi-sensor images and InSAR water level and vegetation monitoring datasets could facilitate wetland change analysis [24,43,44]. Finally, using a different channel of SAR images (e.g., X-, C-, and L-band) in InSAR and classification processes could provide constructive information to characterize coastal wetlands in the Great Lakes and, therefore, reduce misclassification [24,[43][44][45].

Conclusions
Remote sensing has allowed for effective wetland mapping and monitoring applications due to improvements in coverage, accessibility, and processing efficiency. This study utilized the GEE cloud platform and the available Landsat archive to classify wetland and non-wetland classes across the Great Lakes basin over four decades. The resultant classifications achieved an overall accuracy greater than 84% and the change analysis showed that approximately 16% of the study area has undergone landscape change since 1984. The highest increase and decrease in the areas of the classes were observed for Cropland and Forest/Marsh, respectively. Forest mostly transitioned to Fen, Cropland, Marsh, and Swamp, respectively. A notable portion of Marsh was additionally converted into Cropland. These results show that GEE can effectively be applied to monitor wetland change throughout the Great Lakes basin, where continued and improved monitoring will assist in understanding past and present wetland changes for future conservation efforts.