Remote Sensing of Ecosystem Structure—Part 2: Initial Findings of Ecosystem Functioning through Intra- and Inter-Annual Comparisons with Earth Observation Data

: This study examines the response of a cold-regions deltaic wetland ecosystem in northwestern Canada to two separate and differing seasonal wetting cycles. The goal of this paper was to examine the nature of reﬂected electromagnetic energy measured by earth observation (EO) satellites, and to assess whether seasonal wetland hydroperiod and episodic ﬂooding events impact the information retrieved by the Sentinel-2 sensors. The year 2018 represents a year characterized by a large spring freshet and ice-jam ﬂooding, while 2019 represents a year characterized more by summer open-water ﬂooding. We applied the Modiﬁed Normalized Difference Wetness Index (MNDWI) to address the effects of the wetting cycles. The response of the vegetative cover was tracked using the fraction of the absorbed photosynthetically active radiation (fAPAR) and the Leaf Area Index (LAI). All three indices were viewed through the lens of cover classes as derived through a previously published study by the authors. The study provides a framework for designing longer-term studies where multiple intra- and inter-annual hydrological cycles can be accessed via EO data. Future studies will enable the examination of lag times inherent in the response to the various water sources applied to spectral response and incorporate this EO approach into a monitoring framework. These results demonstrate the utility of using indices such as MNDWI, fAPAR and LAI to track vegetation responses in a complex, cold-regions deltaic wetland ecosystem to two separate and differing seasonal wetting cycles. Future work utilizing a longer record, incorporating various combinations of large vs. small snowmelt years, as well as both ﬂooding and non-ﬂooding years, should provide a greater insight into the biological processes across the PAD. As the Sentinel-2 data archive grows, there will be an opportunity to expand this work to include multiple wetting and drying cycles, as well as examine lag times inherent in the response to various hydroclimatic conditions. Integrating the EO-based methodology outlined in this paper into an environmental monitoring framework to support ecosystem assessments, such as for the Wood Buffalo National Park World Heritage Action Plan [46] and Canada-Alberta Oil Sands Monitoring Program [47], will provide an opportunity to assess potential cumulative effects from regional climate change/variability, ﬂow regulation and mining on the downstream PAD region ecosystem.


Introduction
Wetlands are important features in the landscape that provide numerous ecosystem services for people and animals, including storage of water and carbon, protecting and improving water quality, providing essential habitats and access to food, hunting and trapping areas, etc. [1]. These valuable functions are the result of the unique natural characteristics of wetlands; namely, that they form where the landscape is permeated or submerged by water on a permanent or temporary basis for a sufficient length of time to promote aquatic processes. With vegetation adapted to fluctuating water level/depth, wetlands are among the most productive ecosystems in the world, and the only ecosystem designated for conservation by an international convention [2]. Wetlands represent sensitive transition areas between terrestrial and aquatic ecosystems and thus, provide a critical indicator of environmental health [3]. Monitoring of the state of water extent and natural vegetation cover, along with associated processes such as green-up and senescence, provides us with insight as to the stresses being exerted on these sensitive areas as well as their state of health [4].
Earth observation (EO) data have long been applied to the monitoring of various landscapes and ecosystems. Toolboxes commonly applied to EO data contain a myriad of vegetation indices that have been used to assess the state of the vegetation cover [5,6]. Many of these indices are designed to use narrowly focused band-sets (known as narrow band The river-lake-delta system drains Cordilleran and Boreal areas of approximately 600,000 km 2 flowing northwards towards the Slave River. During high water events on the lower Peace River resulting from ice jams and high flows, drainage may reverse southward into the central lakes, and in extreme events, via overland pathways due to small elevational gradients [19,20]. The Peace River has been regulated since 1968 by a number of dams, which control the headwater runoff for the generation of hydroelectricity, classes are generally organized along an elevation and moisture gradient: from forests (balsam poplar, white spruce) at the highest elevations, followed by thickets and savannahs (willows, grasses), marshes and meadows (sedges, grasses, rushes, cattails and reeds), and aquatic communities at the lowest elevations (algae-diatoms, bladderwortduckweed-coontail, pond-lily) [29]. An example of such vegetation cover gradient and wetland channel complexes in the study area of interest within the Athabasca Delta are shown in Figures 2 and 3.

Data and Processing
Sentinel-2 A and B data were obtained for the years 2018 and 2019. From the available image archive, largely cloud-free scenes spanning the growing seasons, and with relatively close calendar date correspondence between the years, were selected for analyses. As such, we acquired images during the open-water season from May, July, August, and September for both years (Table 1). Level 1c processed data, which provided orthorectified TOA (top-of-atmosphere) reflectance data, were employed.   Figure 4, (ii) polygon (green) of areas in Figure 6, and (iii) sampling point design for the current study (yellow polygons). The star indicates the location of picture shown in Figure 2.
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 21 data, and later modified by Du et al. [33] to utilize the longer wavelengths on the ETM+ and Sentinel-2 sensors was selected to describe the land surface wetness. fAPAR and LAI were also chosen as they more closely describe meaningful vegetated conditions [15].

Figure 4.
A spatial representation of the structural classes of the cover for the study area. Reproduced with permission from Peters et al. [18].
LAI is defined as half of the total foliage area per unit ground area [34], and has been used extensively as a proxy for physical and biological plant processes [35,36]. Sellers [37] noted that there was a maximum PAR absorption at LAI of 2 to 3, with absorption saturation >3 resulting in no change with changes in LAI. Sellers [37] also notes that there is some correspondence between area-based physiological processes and multi-spectrally based indices, but that the relationships are attenuated by effects of vegetation species, leaf orientation, LAI, and amount of plant litter. However, a low LAI was interpreted as an indicator of area physiological response to environmental inputs.
VIs used in this work were derived through the Sentinel Application Platform (SNAP) tool as published by the European Space Agency (ESA). The SNAP algorithm, used to derive the indices, was calibrated to rely on Level 1c or TOA reflectance data for the derivation of the indices. The algorithms behind each of the indices are described by Weiss and Baret [9].

Hydrometeorological Data
To evaluate the relative differences in hydrology and meteorology influencing the PAD for the two years of image sequences, flow on the lower Athabasca River at the Embarras Airport hydrometric station (07DD001) entering the AOI were obtained from the Water Survey of Canada [38] and the air temperature and precipitation data from the Meteorological Service of Canada for the nearby climate station at Fort Chipewyan (~30-40 km to the north of the AOI) available through the Government of Alberta portal [39]. A 1 st of April snowpack index was calculated from accumulating the daily total precipitation from 1 November to 31 March [22]. The date of the spring 0 °C isotherm was calculated as the first positive air temperature along a 31-day running means [40]. Potential evapotranspiration was estimated using the Hamon [41] model based on daylength and air temperature that was validated for wetlands in the AOI by Peters et al. [27]. The river-lake-delta system drains Cordilleran and Boreal areas of approximately 600,000 km 2 flowing northwards towards the Slave River. During high water events on the lower Peace River resulting from ice jams and high flows, drainage may reverse southward into the central lakes, and in extreme events, via overland pathways due to small elevational gradients [19,20]. The Peace River has been regulated since 1968 by a number of dams, which control the headwater runoff for the generation of hydroelectricity, leading to important changes to the flow regime near and within the delta, which were partially remediated with the building of weirs on outflow channels [19][20][21]. The Athabasca River mainstem is unregulated, with <5% of the natural flow allocated for agriculture, municipal and industrial water uses and has, thus, retained a near natural hydrological regime [22,23]. Oil sanding mining (open pit and in situ) is occurring in tributaries of the Lower Athabasca River, just upstream of the PAD, where landscape alterations and small amounts of water are abstracted from the mainstem for processing of bitumen [23,24].
The PAD contains >1000 inland wetland and lake basins with varying degrees of connectivity to the main flow system: open drainage, restricted and isolated basins [25,26]. The latter two classes are perched above and separated from the main flow system. Other than areas with significant local runoff input (e.g., Shield areas in northeast area of the PAD), restricted and isolated basins require periodic replenishment of water via overland flow routes activated during spring breakup ice jams or summer high flows to offset loses of water in between inundation events mainly through potential evapotranspiration. Potential evapotranspiration rates are, on average, greater (~10 cm) than local precipitation [27,28]. The generally prevailing semi-arid climate dominated by evapotranspiration rates, therefore, strongly influences surface water extents and depths in this area, with periodic flooding and water drawdown cycle resulting in a spatiotemporally variable wetland landscape, but highly productive ecosystem [29].
The PAD has been recognized as a primary habitat for a large range of insects, birds and mammals [29][30][31]. The vegetation cover, as reported on by Timoney [30], is affected by topography, frequency of flooding and depth to the groundwater table. Vegetation classes are generally organized along an elevation and moisture gradient: from forests (balsam poplar, white spruce) at the highest elevations, followed by thickets and savannahs (willows, grasses), Remote Sens. 2021, 13, 3219 6 of 21 marshes and meadows (sedges, grasses, rushes, cattails and reeds), and aquatic communities at the lowest elevations (algae-diatoms, bladderwort-duckweed-coontail, pond-lily) [29]. An example of such vegetation cover gradient and wetland channel complexes in the study area of interest within the Athabasca Delta are shown in Figures 2 and 3.

Data and Processing
Sentinel-2 A and B data were obtained for the years 2018 and 2019. From the available image archive, largely cloud-free scenes spanning the growing seasons, and with relatively close calendar date correspondence between the years, were selected for analyses. As such, we acquired images during the open-water season from May, July, August, and September for both years (Table 1). Level 1c processed data, which provided orthorectified TOA (top-of-atmosphere) reflectance data, were employed.

Ecosystem Structure
The study area was originally classified based on an integration of Sentinel-2 imagery and airborne LiDAR data into reflectance/structurally based cover classes. A detailed description of the process by which the cover classes were generated is provided in Peters et al. [18]. A visual, spatial summary of the segmentation of the landscape is provided in Figure 4, with the cluster numbering sequence retained in this study for the description of the behavior of the time series-based VI.

Spectral Indices
As a precursor to our analysis of the consistency of the annual vegetative response, the degree and persistence of wetting conditions throughout the 2018 and 2019 growing seasons were examined. The MNDWI originally proposed by McFeetters [32] for LANDSAT MSS data, and later modified by Du et al. [33] to utilize the longer wavelengths on the ETM+ and Sentinel-2 sensors was selected to describe the land surface wetness. fAPAR and LAI were also chosen as they more closely describe meaningful vegetated conditions [15].
LAI is defined as half of the total foliage area per unit ground area [34], and has been used extensively as a proxy for physical and biological plant processes [35,36]. Sellers [37] noted that there was a maximum PAR absorption at LAI of 2 to 3, with absorption saturation >3 resulting in no change with changes in LAI. Sellers [37] also notes that there is some correspondence between area-based physiological processes and multi-spectrally based indices, but that the relationships are attenuated by effects of vegetation species, leaf orientation, LAI, and amount of plant litter. However, a low LAI was interpreted as an indicator of area physiological response to environmental inputs.
VIs used in this work were derived through the Sentinel Application Platform (SNAP) tool as published by the European Space Agency (ESA). The SNAP algorithm, used to derive the indices, was calibrated to rely on Level 1c or TOA reflectance data for the derivation of the indices. The algorithms behind each of the indices are described by Weiss and Baret [9].

Hydrometeorological Data
To evaluate the relative differences in hydrology and meteorology influencing the PAD for the two years of image sequences, flow on the lower Athabasca River at the Embarras Airport hydrometric station (07DD001) entering the AOI were obtained from the Water Survey of Canada [38] and the air temperature and precipitation data from the Meteorological Service of Canada for the nearby climate station at Fort Chipewyan (~30-40 km to the north of the AOI) available through the Government of Alberta portal [39]. A 1st of April snowpack index was calculated from accumulating the daily total precipitation from 1 November to 31 March [22]. The date of the spring 0 • C isotherm was calculated as the first positive air temperature along a 31-day running means [40]. Potential evapotranspiration was estimated using the Hamon [41] model based on daylength and air temperature that was validated for wetlands in the AOI by Peters et al. [27].

Comparison of Wetting Conditions for the 2018 and 2019 Growing Seasons
Compared to the 30-year mean climate, both study years 2018 and 2019 were slightly cooler than normal, received less than normal annual precipitation, and experienced above and below normal rates of total potential evapotranspiration, respectively ( Table 2). Noteworthy to our analysis timeline was the spring ice-jam flood event of 2018 (later date of snowmelt as indicated by 0 • C isotherm, greater 1st of April snowpack index, and a larger spring peak flow magnitude coupled with mechanical ice breakup) around the lower Athabasca River, Embarras River and Mamawi Creek region that inundated adjacent lakes and wetlands in the delta floodplain covered by our AOI [42] ( Figure 5). Furthermore, both summers experienced open-water flooding/connectivity events from high flows (i.e., >1750 m 3 s −1 ) entering the Athabasca Delta as observed in several wetlands along the Mamawi Creek (unpublished data, Peters), with 2019 experiencing large and multiple summer overbanking events. Although the year 2019 had an earlier start to spring 0 • C isotherm, 2018 experienced considerably more potential evapotranspiration and less precipitation, leading to greater climatological net water loss during the growing season than in 2019 (i.e., precipitation minus potential evapotranspiration).  Figure 6 illustrates the spatial and temporal variations in the wetness conditions for the two years as determined by the MNDWI. As is evident in the two annual sequences of imagery, the 2018 data ( Figure 6a-d) indicate that the May through August periods were characterized by wetter conditions than in the September scene. The levee and floodplain areas bordering the Athabasca River and Embarras River channels in the southern portion of the study area indicate extensive wetness. The intensity of this wetness decreases in July, increasing in the August and decreases in the September 2018 MNDWI image in the areas bordering the Athabasca River in the south. This is especially apparent in the early scenes, while the August scene indicates that there was a widespread increase in the MNDWI. If we relate the local temporal hydroclimatic patterns in Figure 5a to these MNDWI images, we can see that there were inundation events and several periods of rainfall prior to the acquisition dates for the three first scenes. There was a relatively dry period and declining flows leading up to the September acquisition, which is also reflected in the MNDWI output for that date.   The 2019 MNDWI imagery (Figure 6e-h) provides us with a different scenario, with an initially drier spring characterized with relatively lower spring flows (i.e., no ice-jam flooding), followed by a wetting period through August and then drying into the start of Autumn. These patterns are consistent with the hydrometeorological patterns presented in Figure 5b above. For the 2019 growing season, there was a period of precipitation which started prior to the July acquisition date, lasting throughout the summer until after the August imaging date. Moisture was also diverted into areas of the delta landscape as a result of several high flow events. There was an extended drying period prior to the September image.
For a more detailed examination of the inter-and intra-year changes, we generated a set of 4000 randomly spaced points on which to base our sample. The distribution of these points is from south to north (see Figure 3 for polygons of points sampled), which corresponds to the observable variations in moisture content, especially as related to the 2018 flooding. These points were chosen because the areas bounded by the polygons: (i) encompassed the range of cover classes present in the study area; (ii) minimized the occurrence of standing water in the form of shallow ponds and lakes; and (iii) known to have experienced a range of flood timing and flood types. These same points were subsequently used below to address the changing vegetation indices. The intra-year comparison of the wetting patterns provided in Figure 7 shows that for 2018, there were inconsistencies between the early images. In comparing the July-August data, there was a general wetting of the area, while the August-September comparison was more organized, although the area was slightly wetter in the August scene. For 2019, the sample point comparison supports our more qualitative interpretation of Figure 6. Note that a 1:1 line was added to Figures 7 and 8 to facilitate interpretation; in the case of no difference between the two periods, the points would lie on this line. image in the areas bordering the Athabasca River in the south. This is especially apparent in the early scenes, while the August scene indicates that there was a widespread increase in the MNDWI. If we relate the local temporal hydroclimatic patterns in Figure 5a to these MNDWI images, we can see that there were inundation events and several periods of rainfall prior to the acquisition dates for the three first scenes. There was a relatively dry period and declining flows leading up to the September acquisition, which is also reflected in the MNDWI output for that date.  The inter-year comparison presented in Figure 8 again supports our previous interpretation of Figure 6, with 2018 being considerably wetter in May, August and September than for the same periods in 2019 as indicated by MNDWI. Of the months examined, July appeared to be the more consistent between the two study years, but that could be attributed to the hydrometeorology of these specific years and may not extend to other years. Overall, wetness in the AOI of the Athabasca Delta was highly dynamic in response to various water inputs and outputs.
sponds to the observable variations in moisture content, especially as related to the 2018 flooding. These points were chosen because the areas bounded by the polygons: (i) encompassed the range of cover classes present in the study area; (ii) minimized the occurrence of standing water in the form of shallow ponds and lakes; and (iii) known to have experienced a range of flood timing and flood types. These same points were subsequently used below to address the changing vegetation indices. The intra-year comparison of the wetting patterns provided in Figure 7 shows that for 2018, there were  While meteorological-wetness associations would account for a general wetting or drying of the study area, they do not account for the spatial wetting/drying patterns. If we examine the May 2018 MNDWI image (Figure 6a-d), we see that the wetness was concentrated around the Athabasca and Embarras Rivers in the southern portions of the AOI, which is an important inflow zone of runoff generated from >150,000 km 2 , while the farther inland and lower elevation areas in the more northern deltaic portions of the study area remained dry relative to the higher areas closer to these rivers. This pattern is not as apparent in the time sequence of images for 2019 (Figure 6e-h). This is consistent with the 2018 flooding event, which was a result of ice jams temporarily impeding incoming spring freshet flow and generating extremely high waters that flowed over banks during that event. The year 2019, which provides us with a very different spatial wetting and drying pattern is more consistent with what we would expect with no spring ice jam, and summer rainfall-related processes enhanced by localized high flows overbanking channels into adjacent deltaic landscape, followed by drainage and drawdown.
In addition to the general wetting and drying pattern within the terrestrial portion of the delta, there were changes to the number and surface water extent of wetlands and small lakes scattered throughout the AOI. In examining Figure 6, we can see that in 2018, the small lakes and open water wetlands in the top half of the study area were larger and more numerous. The 2019 pattern is different in that the water bodies were less numerous in the early (May) image in comparison to July, which is particularly evident in the areas between the Athabasca and Embarras Rivers. MNDWI values surrounding these water bodies are lower than seen in the later data, indicating that not only has there been draining subsequent to inundations, but also a general drying of the landscape throughout the course of the growing season via evapotranspiration. a general wetting of the area, while the August-September comparison was more organized, although the area was slightly wetter in the August scene. For 2019, the sample point comparison supports our more qualitative interpretation of Figure 6. Note that a 1:1 line was added to Figures 7 and 8 to facilitate interpretation; in the case of no difference between the two periods, the points would lie on this line.
The inter-year comparison presented in Figure 8 again supports our previous interpretation of Figure 6, with 2018 being considerably wetter in May, August and September than for the same periods in 2019 as indicated by MNDWI. Of the months examined, July appeared to be the more consistent between the two study years, but that could be attributed to the hydrometeorology of these specific years and may not extend to other years. Overall, wetness in the AOI of the Athabasca Delta was highly dynamic in response to various water inputs and outputs. While meteorological-wetness associations would account for a general wetting or drying of the study area, they do not account for the spatial wetting/drying patterns. If we examine the May 2018 MNDWI image (Figure 6a-d), we see that the wetness was concentrated around the Athabasca and Embarras Rivers in the southern portions of the AOI, which is an important inflow zone of runoff generated from >150,000 km 2 , while the farther inland and lower elevation areas in the more northern deltaic portions of the study area remained dry relative to the higher areas closer to these rivers. This pattern is not as apparent in the time sequence of images for 2019 (Figure 6e-h). This is consistent with the 2018 flooding event, which was a result of ice jams temporarily impeding incoming spring freshet flow and generating extremely high waters that flowed over banks during that event. The year 2019, which provides us with a very different spatial wetting and drying pattern is more consistent with what we would expect with no spring ice jam, and summer rainfall-related processes enhanced by localized high flows overbanking channels into adjacent deltaic landscape, followed by drainage and drawdown.

Response Comparison of Vegetation-Related Spectral Indices for 2018 and 2019
Capitalizing on the spatiotemporally changing hydrometeorological conditions in the study area, and the interdependency of the vegetation and the hydroperiod and moisture within wetland/lake basins and surrounding landscape, we examined the response of the vegetation indices during the 2018 and 2019 snow/ice-free seasons. To make this assessment, the same random point set described above was used to address the indices (LAI, fAPAR, MNDWI) from both years corresponding to the cover class cluster number found in Figure 4.
We performed a test of the separability of the clusters using a Kruskal-Wallis test. The results of these statistical tests on an intra-and inter-year basis are presented in Appendix A Tables A1-A4. For the sake of clarity, the individual clusters were grouped into their interpreted cover classes (see Figure 4), and the examination was restricted to the non-water classes. This was done due to the transient nature of the water bodies over time and that LAI and fAPAR for water bodies may not be the most appropriate metrics to use. Kruskal-Wallis tests applied to the pooled samples indicated that there are no statistical differences between the datasets for the same periods for the two years (Table A2). The exception being for the fAPAR between May and July 2018, with the remaining two periods for 2018 showing no significant differences (p > 0.05). Consistent with observed fAPAR values in 2019, the three comparison periods for LAI returned no significant differences in the observed output. Table A4 summarizes the examination of the changes in the within-year cluster position and in all cases, the LAI and fAPAR index clusters differed seasonally.
A summary of the LAI and fAPAR for both the 2018 and 2019 datasets are presented in Figures 9-12. The behavior of the fAPAR and LAI during 2018 for the coniferous class is relatively consistent between periods (Figure 9). The early (May) fAPAR data cluster is slightly lower than that seen for the end of growing seasons (September) data. The July and August fAPAR clusters are at the same levels, even though the MNDWI levels for August are lower (drier). A similar pattern is observed with the 2018 LAI data with MNDWI. For the May 2019 fAPAR cluster, the levels are lower than we see for the other months. The fAPAR data for the three remaining 2019 acquisition dates are relatively consistent, although the MNDWI data ranges are not. The LAI data suggest that there was an increase in the overall LAI data from May through the summer acquisitions. However, the September cluster LAI levels decreased and are somewhat consistent with those seen in the May imagery. The differences between the May and September are primarily in the MNDWI levels. When we examine the relationships between fAPAR and LAI, we see a consistency between the May, July and August dates, with higher fAPAR values relative to LAI in September. These patterns can also be seen in the 2019 datasets, although there is less overlap between the May and July/August data.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 21 but also maintaining a similar level into September. LAI patterns were similar again, with the May imagery being lower than the other three periods. The Emergent vegetation ( Figure 12) behaved differently than for the other three classes observed in that the 2018 early (May) fAPAR data are consistent with the observations for July and September, and that there is little change in this index for these three observation dates. The LAI data indicate less overlap between the early/late and mid-season values. For the 2019 fAPAR data, the pattern observed is more consistent with the other cover classes; however, the clustering is not as tight in this class as in the other three.  The Deciduous and Herbaceous/Meadow class patterns (Figures 10 and 11) are similar to the those observed for the Coniferous class for the 2018 fAPAR data. There is a gradual increase in fAPAR from May through July and August, with a decrease into September. A similar pattern can be seen for the LAI values. In 2019, the fAPAR metrics were lower in the May image than we see for 2018, rising to equivalent levels in July and August but also maintaining a similar level into September. LAI patterns were similar again, with the May imagery being lower than the other three periods.
The Emergent vegetation ( Figure 12) behaved differently than for the other three classes observed in that the 2018 early (May) fAPAR data are consistent with the observations for July and September, and that there is little change in this index for these three observation dates. The LAI data indicate less overlap between the early/late and midseason values. For the 2019 fAPAR data, the pattern observed is more consistent with the other cover classes; however, the clustering is not as tight in this class as in the other three.     All three classes suggest a much lower level of development/green-up than is the case for the Coniferous class in 2019. The July and August data both yield similar patterns for both dates. The September data indicate slightly higher fAPAR levels in 2018, but also less organized fAPAR levels at lower LAI. This may suggest the initial onset of senescence for this later near Fall date, when mean daily air temperatures dropped below 5 • C and daylength shortened to <12 h.

Discussion
The wetting patterns observed through the MNDWI spatiotemporal sequences presented in Figure 6 indicate that water dispersal in a complex, cold-regions wetland-lake landscape can be tracked through the use of EO data such as those acquired by the Sentinel-2 platform. The observed patterns in the 2018 dataset clearly indicate that the spring wetting radiated out from the major river and distributary channels that feed into the Athabasca Delta sector of the PAD. It is interesting to note that the areas to the north of our study area, which are generally at lower elevations and in closer proximity to the large central lakes (Lake Athabasca, Mamawi Lake, Lake Claire), were distant from and subsequently not as affected by ice-jam flooding, and thus, seen as drier than the southern areas more proximal to the river channels. It is also notable in the 2018 July and August data that the image acquisitions follow periods of rainfall and high flow events which would account for the increasing wetness, and that there was a drying pattern in the September image date consistent with the period of little rainfall immediately preceding the acquisition.
The wetting signal for 2019 was very different given that there was a relatively small spring runoff and no ice jams to divert water onto the landscape early in the growing season. The May image reveals a relatively dry study area that was influenced mostly by a local, below normal snowmelt amount. Wetting occurred throughout July and August consistent with a wetter mid-summer from higher rainfall and high flow-driven inundations of areas proximal to the river channels (e.g., Mamawi Creek). However, these areas were not as spatially pronounced as those affected by an ice-induced flood mechanism, which was capable of shunting water farther inland and into more elevated floodplain areas [19]. Similar to 2018, a general drying period prior to the September image date was seen in 2019.
The observed patterns in Figures 9-12 are generally consistent with those observed by Kozlov et al. [43] who note that for large boreal wetlands in Russia, there is a general increase in fAPAR during the year peaking at around June. The wetter "marsh" vegetation, they note, has fAPAR peaks in the order of 0.8, while the woodland values peak at fAPAR 0.5-0.6. In our study, we see somewhat higher fAPAR levels with the coniferous samples for 2018, but a lower peak in 2019. Deciduous stands are consistent between the two years. Herbaceous samples peak at 0.8 in both 2018 and 2019, which is in line with the observations reported by Kozlov et al. [43]. What is particularly interesting to observe in Figures 9-12 is the consistent cyclical nature of the point distributions patterns. The exception is in the case of the early season (May) observations, which appear to be more consistent with variations in the hydroclimate and flooding data. The patterns of observations from the later parts of the growing season may speak more to the influence of the long periods of solar insolation experienced in northern regions, which may dominate the growth cycle rather than water being the main growth-influencing variable (see expansion of this point below).
The relationship between fAPAR and LAI as described in Figures 13-16 for the woodland samples (coniferous and deciduous) is typically linear for lower LAI values, but when considering the entire range of LAI, values tend towards a nonlinear relationship towards the upper end of the LAI range (LAI > 2.0). This can be interpreted as a decreasing sensitivity of satellite-derived LAI to describe changes in fAPAR values during the midsummer period when the foliage is at its peak. This is consistent with observations from others that have examined the behavior of vegetation indices as they relate to vegetation structure [44,45]. An interesting observation from the data point in Figures 13-16 is that the relationship between LAI and fAPAR is highly correlated. These two indices were chosen as they were thought to represent the canopy reaction to environmental inputs at different time scales (fAPAR more immediate, LAI more delayed). The take-away from this observation is that the LAI-based index could be used as a surrogate for fAPAR and foliar functioning. This is an important detail if we transition to a LANDSAT-based archive, which contains a longer time series, but with fewer spectral bands on which to derive a more limited set of vegetation indices.
As discussed above, the daily hydrometeorological data summarized in Figure 5 provide some insight as to the key environmental variables and their potential impact on the observed EO data. The MNDWI data for May 2018 map a general wetter surface, the cause of which may be related to the snowmelt/rainfall immediately preceding the May acquisition. However, this does not explain the spatial distribution of the MNDWI values observed in Figure 6. The wetness patterns for 2018, especially in the early image, radiate out from the larger river channels from south to north, which is more reflective of the spring flooding episode, rather than weather-related factors. While the MNDWI values for the 2019 map are spatially more homogeneous and indicate a drier overall surface over the growing season, there does not seem to be as dramatic a correlation with the wetter rainfall and peak flow period recorded in the weeks prior to the August image date. The sampled fAPAR and LAI values for those two image dates also indicate a consistency between the mid-summer sample dates for the two years. It is predominantly the May samples that provide a major difference between the two study years. The 2018 data suggest a higher level of photosynthetic activity and leaf area in the May data than observed for the same month in 2019. This is in spite of the fact that the spring melt commenced earlier in 2019 than it did in 2018. However, the May mean air temperature in 2018 was 4 • C higher than in 2019 (10.8 • C vs. 6.2 • C, respectively), promoting greater vegetation growth. The July and August observations suggest that the LAI and fAPAR measurements normalize during the mid-summer period and that the supposed effects of the spring flooding are mitigated, unless the occurrence of subsequent inundation events.

Conclusions
While this work examined one extensive spring ice-jam flooding and one less extensive open-water flooding year, there was an apparent early biological response to the former.
These results demonstrate the utility of using indices such as MNDWI, fAPAR and LAI to track vegetation responses in a complex, cold-regions deltaic wetland ecosystem to two separate and differing seasonal wetting cycles. Future work utilizing a longer record, incorporating various combinations of large vs. small snowmelt years, as well as both flooding and non-flooding years, should provide a greater insight into the biological processes across the PAD. As the Sentinel-2 data archive grows, there will be an opportunity to expand this work to include multiple wetting and drying cycles, as well as examine lag times inherent in the response to various hydroclimatic conditions. Integrating the EObased methodology outlined in this paper into an environmental monitoring framework to support ecosystem assessments, such as for the Wood Buffalo National Park World Heritage Action Plan [46] and Canada-Alberta Oil Sands Monitoring Program [47], will provide an opportunity to assess potential cumulative effects from regional climate change/variability, flow regulation and mining on the downstream PAD region ecosystem.   Table A2. Kruskal-Wallis test results for fAPAR between adjacent image acquisitions dates. The presence of *** denotes that the null hypothesis was rejected and that the samples are significantly different at p = 0.05 level.  Table A3. Kruskal-Wallis test results for LAI between adjacent image acquisitions dates. The presence of *** denotes that the null hypothesis was rejected and that the samples are significantly different at p = 0.05 level.  Table A4. Kruskal-Wallis test results for comparison of separation metrics (LAI and fAPAR) based on vegetation cover class/clusters through time for LAI between adjacent image acquisitions dates. The presence of *** denotes that the null hypothesis was rejected and that the samples are significantly different at p = 0.05 level.