Wetland Hydroperiod Change Along the Upper Columbia River Floodplain, Canada, 1984 to 2019

: Increasing air temperatures and changing hydrological conditions in the mountainous Kootenay Region of British Columbia, Canada are expected to a ﬀ ect ﬂoodplain wetland extent and function along the Columbia River. The objective of this study was to determine the seasonally inundated hydroperiod for a ﬂoodplain section (28.66 km 2 ) of the Upper Columbia River wetlands complex using time series satellite image observations and binary open water mask extraction. A mid pixel resolution (30 m) optical satellite image time series of 61 clear sky scenes from the Landsat Thematic Mapper


Remote Sensing of Water Extent and Hydroperiod
Mapping changing open water extent or hydroperiod over entire seasons or periods of years requires temporally frequent and spatially explicit observations that can only be practically achieved through satellite remote sensing [19,[48][49][50][51]. NASA's Landsat mission has the longest almost continuous record of high-resolution (<100 m pixel resolution) optical satellite image data stretching back to 1972 with the Multi-Spectral Scanner (MSS) sensor, which was operational until 1998 (26 years). However, the higher resolution Thematic Mapper (TM) and now Optical Land Imager (OLI) sensors provide a still active and therefore longer record of observation from 1984 to the present day (>36 years). These long time series lend themselves to evaluating changes in water body extent and hydroperiod mapping. Several studies have demonstrated optical water mask extraction using spectral index or normalized ratio approaches [52][53][54][55]. For example, the difference between the Normalized Difference Vegetation Index (NDVI) and the Modified Normalized Water Index (MNDWI) has been found effective for water mask extraction as it effectively removes interference, enhances the differentiation between land and water and is well suited for areas of complex topography. Due to the orbital repeat of 16 days and image occlusions during periods of cloud cover, however, Landsat images are unsuited to hydroperiod mapping during a single year.
Due to its active sensing capability at microwave wavelengths, synthetic aperture radar (SAR) imaging is ideal for mapping surface water extent [56][57][58][59] at any time of day or night and through cloud cover. The specular reflectance of water causes little backscatter to the sensor [60], allowing water boundary extraction algorithms to be developed using a variety of approaches [59,[61][62][63]. Airborne lidar (light detection and ranging) is also effective at mapping wetland water bodies [64] and features [65,66] at very high resolutions. Data fusion using multiple remote sensing sources has resulted in improved discrimination of terrain, vegetation and water body information to enhance wetland classification and hydrologic connectivity [65,67]. For example, many studies have shown that lidar terrain information can play a key role in the quality control of water masks derived from other sensing technologies by identifying areas that cannot be inundated by virtue of their elevation or terrain morphology [49][50][51]68].

Objectives
The underlying premise for this study is that if changes in floodplain wetland hydroperiod can be discerned and related to historical changes in regional hydroclimatic conditions, then potential future climate change impacts on wetland and water body permanence can be inferred. As part of a broader wetland vulnerability assessment framework, therefore, the primary objective is to characterize the multidecadal wetland complex hydroperiod for a section of the Upper Columbia River floodplain from the Landsat TM and OLI archive from 1984 to 2019. This 36-year hydroperiod can be divided into an early and late period so that any change in the seasonal extent and permanence of open water wetlands can be evaluated. Underlying drivers of any observed shift in hydroperiod behaviour will be investigated by extracting nearby meteorological norms for the two periods and testing for a relationship between downstream river flow and individual open water mask extents.
A sub-objective is to evaluate the potential for extending the hydroperiod backwards using lower spatial resolution archive Landsat MSS data and forwards through continued monitoring using high temporal and spatial resolution SAR imagery, such as is becoming available through the Radarsat Constellation Mission (RCM) [69]. The 4-day repeat cycle of RCM will provide high temporal resolution time series which will facilitate hydroperiod determination.

Study Area
The case study wetland complex is found along a~17 km stretch of the Upper Columbia River floodplain between Brisco and Spillimacheen (50 •  to downstream end ranging from 789 to 785 m a.s.l. (hydraulic gradient of 0.02%). Nestled between high mountain ranges (up to~3500 m a.s.l.), including the Rocky Mountains to the east and the Columbia and Purcell Mountains in the west, the deep post-glacial valley is fed by multiple tributaries [28]. Surrounded by forests, agriculture, glaciers and park land, the Upper Columbia River valley hosts many wetland complexes which are often hydrologically disconnected from the main channel and considered potentially vulnerable to a warming climate [27,28].
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 20 downstream end ranging from 789 to 785 m a.s.l. (hydraulic gradient of 0.02%). Nestled between high mountain ranges (up to ~3500 m a.s.l.), including the Rocky Mountains to the east and the Columbia and Purcell Mountains in the west, the deep post-glacial valley is fed by multiple tributaries [28]. Surrounded by forests, agriculture, glaciers and park land, the Upper Columbia River valley hosts many wetland complexes which are often hydrologically disconnected from the main channel and considered potentially vulnerable to a warming climate [27,28]. Summer daytime high air temperatures range from 14 to 19 °C, with winter daytime high temperatures of 4-8 °C. Lower elevations receive less precipitation (325-471 mm/year) than higher elevations (991-1531 mm/year), where the majority of precipitation occurs during winter months [28]. The Columbia River is gauged continuously by the Water Survey of Canada at Nicholson, ~65 km downstream from the case study site, where the upstream drainage area is ~6660 km 2 . During the 1984 to 2019 study period, annual peak flows tended to occur during June and July (Figure 2), ranging from 282 (2001) to 743 m 3 s −1 (2002). Annual base flow conditions are relatively stable between 25 and 35 m 3 s −1 and tend to be present from late October through to early April. Floodplain wetlands connected to the mainstem of the Columbia River experience water level and extent fluctuations as river flow rises and falls, while more isolated (less connected) wetland water bodies experience level and extent changes as a result of overbank flooding, slow drainage and loss to evapotranspiration [28,43]. Summer daytime high air temperatures range from 14 to 19 • C, with winter daytime high temperatures of 4-8 • C. Lower elevations receive less precipitation (325-471 mm/year) than higher elevations (991-1531 mm/year), where the majority of precipitation occurs during winter months [28]. The Columbia River is gauged continuously by the Water Survey of Canada at Nicholson,~65 km downstream from the case study site, where the upstream drainage area is~6660 km 2 . During the 1984 to 2019 study period, annual peak flows tended to occur during June and July (Figure 2), ranging from 282 (2001) to 743 m 3 s −1 (2002). Annual base flow conditions are relatively stable between 25 and 35 m 3 s −1 and tend to be present from late October through to early April. Floodplain wetlands connected to the mainstem of the Columbia River experience water level and extent fluctuations as river flow rises and falls, while more isolated (less connected) wetland water bodies experience level and extent changes as a result of overbank flooding, slow drainage and loss to evapotranspiration [28,43].

Data Sources
Daily discharge data ( Figure 2) for the Columbia River at Nicholson are publicly available and were downloaded from the Environment and Climate Change Canada (ECCC) online data portal for the years 1984-2019. Similarly, air temperature and precipitation data for Cranbrook BC (~100 km south of the study site) are publicly available from the ECCC online portal and were mined to provide an index of changing climate conditions within the study area. While some distance away, data from Cranbrook were used in this study due to the uninterrupted data availability throughout the 1984 to 2019 period. A 1 m lidar digital elevation model (DEM) was obtained over the Columbia River floodplain by the BC Province and made available to this study for orthorectification of the SAR images, to define the spatial extent of the floodplain and thus constrain the extent of the floodplain water masks. Floodplain extent was manually digitized from the DEM using the elevation contours and break of slope between the flat valley bottom and the surrounding mountain slopes. Given that the DEM resolution is much higher than any of the image data sources used, any errors in digitization will have a negligible impact on water mask quality control.
The publicly accessible USGS EarthExplorer data archive was mined for all cloud-free Landsat TM and OLI scenes over the study site (Path 43, Row 25) during the hydrologically active season from April to October 1984 to 2019. In total, 61 scenes (49 TM and 12 OLI) met the criteria and were downloaded at 30 m pixel resolution as a Collection Level-1 GeoTIFF data product. Of the 35 years sampled, only five years displayed no suitable scenes: 1993, 1999, 2010, 2012 and 2013. For the years of available data, there was an average of two images a year, with a maximum of five in 1994, followed by four in years 2003, 2007, 2016 and 2017. Single scenes were captured during 1984, 1985, 1989, 1990-1992, 1995-1996, 2000, 2006, 2011, 2014. The frequency of suitable images by month from April to October was 3, 3, 4, 16, 20, 14, 1, respectively, with the highest sample frequency during mid to late summer flow conditions and the lowest frequencies close to the early and late season baseflow periods. An additional 12 Landsat MSS scenes were acquired for 1992-1995, where the imaging frequency was, respectively 4, 1, 5, 2 for each of the years.
Recent Radarsat 2 images were obtained from the Canada Centre for Mapping and Earth Observation (CCMEO) of Natural Resources Canada. Of all Radarsat products available, only single HH polarization (Horizontal transmit/Horizontal receive) descending orbit images in Single Look Complex (SLC) format were used, as the high relief terrain surrounding the study area caused shadowing in ascending track scenes and excessive noise in scenes with a vertical polarization component. Only 8 suitable images were available for 2015 to 2019, with, respectively, 1, 4, 2, 0, 1 images for each of the years. Radarsat 2 extra-fine products have a 3 m × 3 m resolution. Of the active season months, 6 of the 7 were sampled at least once, with 2 images each for April and July and no image for September.

Data Sources
Daily discharge data ( Figure 2) for the Columbia River at Nicholson are publicly available and were downloaded from the Environment and Climate Change Canada (ECCC) online data portal for the years 1984-2019. Similarly, air temperature and precipitation data for Cranbrook BC (~100 km south of the study site) are publicly available from the ECCC online portal and were mined to provide an index of changing climate conditions within the study area. While some distance away, data from Cranbrook were used in this study due to the uninterrupted data availability throughout the 1984 to 2019 period. A 1 m lidar digital elevation model (DEM) was obtained over the Columbia River floodplain by the BC Province and made available to this study for orthorectification of the SAR images, to define the spatial extent of the floodplain and thus constrain the extent of the floodplain water masks. Floodplain extent was manually digitized from the DEM using the elevation contours and break of slope between the flat valley bottom and the surrounding mountain slopes. Given that the DEM resolution is much higher than any of the image data sources used, any errors in digitization will have a negligible impact on water mask quality control.
The publicly accessible USGS EarthExplorer data archive was mined for all cloud-free Landsat TM and OLI scenes over the study site (Path 43, Row 25) during the hydrologically active season from April to October 1984 to 2019. In total, 61 scenes (49 TM and 12 OLI) met the criteria and were downloaded at 30 m pixel resolution as a Collection Level-1 GeoTIFF data product. Of the 35 years sampled, only five years displayed no suitable scenes: 1993, 1999, 2010, 2012 and 2013. For the years of available data, there was an average of two images a year, with a maximum of five in 1994, followed by four in years 2003, 2007, 2016 and 2017. Single scenes were captured during 1984, 1985, 1989, 1990-1992, 1995-1996, 2000, 2006, 2011, 2014. The frequency of suitable images by month from April to October was 3, 3, 4, 16, 20, 14, 1, respectively, with the highest sample frequency during mid to late summer flow conditions and the lowest frequencies close to the early and late season baseflow periods. An additional 12 Landsat MSS scenes were acquired for 1992-1995, where the imaging frequency was, respectively 4, 1, 5, 2 for each of the years.
Recent Radarsat 2 images were obtained from the Canada Centre for Mapping and Earth Observation (CCMEO) of Natural Resources Canada. Of all Radarsat products available, only single HH polarization (Horizontal transmit/Horizontal receive) descending orbit images in Single Look Complex (SLC) format were used, as the high relief terrain surrounding the study area caused shadowing in ascending track scenes and excessive noise in scenes with a vertical polarization component. Only 8 suitable images were available for 2015 to 2019, with, respectively, 1, 4, 2, 0, 1 images for each of the years. Radarsat 2 extra-fine products have a 3 m × 3 m resolution. Of the active Remote Sens. 2020, 12, 4084 6 of 20 season months, 6 of the 7 were sampled at least once, with 2 images each for April and July and no image for September.

Image-Based Water Mask Extraction
The workflow for optical and SAR water mask and hydroperiod construction is summarized in Figure 3. For the Landsat images, an index-based threshold water mask was generated using the NDVI and MNDWI indices [55]. Images were radiometrically and atmospherically corrected using the Radiometric Calibration and FLAASH tools in ENVI. Following the Landsat TM band number convention, NDVI layers were created in ENVI using the near infrared (NIR) (4) and red (3) bands in the band math tool as follows: Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 20

Hydroperiod Frequency Analysis
For comparative hydroperiod analysis, the Landsat TM/OLI water masks were divided into early and late time periods (T1 and T2) at the annual increment that produced the closest to an equal number of image mask samples, such that 1984 to 2002 (T1: 18 years, n = 30) was compared to 2003 to The MNDWI layers were created using the green (2) and the short-wave infrared (SWIR) (5) bands using the following equation: (Note, for Landsat OLI, the NDVI calculation uses band 5 and band 4, the MNDWI band 3 and 6, and for Landsat MSS, change the band order to green (1), red (2), NIR (3) and SWIR I (4)). The MNDWI layer was then subtracted from NDVI to provide an image with values ranging between +2 and −2.
The index difference images were then passed through a majority filter with a kernel set to the nearest four orthogonal pixels to reduce pixel-level noise. Chosen water threshold values ranged between 0.15 and 0.19 and were trained for each image over known permanent water bodies throughout the floodplain and manually applied to produce a binary mask where a value of 1 signified water and 0 everything else. The mask was then clipped to the floodplain extent.
SAR-derived surface water masks were created following the intensity (dB) thresholding routine and filter settings of Montgomery et al. [51]. Radarsat 2 images were orthorectified using PCI Ortho Engine (PCI Geomatica) and then imported into PCI Model Builder, where the PSIOInterp tool converted pixel intensity to decibel values. These images were then run through a 3 × 3 average (mean) filter to reduce speckle and a 5 × 5 FGAMMA (gamma maximum-a-posteriori (MAP)) adaptive filter to preserve edges, crucial for water extraction. Both outputs were then manually inspected, and a threshold was trained over known permanent water bodies for each individual image. The range of pixel thresholds was ±14-17 dB, with negative values for the FGAMMA threshold and positive values for the average (mean) filter. Water values were converted to 999 for both images and the two images were then added together. Outputted pixel values ranged from 0 to 1998, with 1998 values identified as water by both filters. This image was then added to the original FGAMMA-filtered image and water pixels assigned a value of 1 and everything else 0. The binary mask was passed through a 3 × 3 majority (or modal) filter to reduce noise and a 9 × 9 FGAMMA filter to preserve edges. As with the optical processing steps, the final water masks were then clipped to the floodplain extent.

Hydroperiod Frequency Analysis
For comparative hydroperiod analysis, the Landsat TM/OLI water masks were divided into early and late time periods (T1 and T2) at the annual increment that produced the closest to an equal number of image mask samples, such that 1984 to 2002 (T1: 18 years, n = 30) was compared to 2003 to 2019 (T2: 16 years, n = 31). The active season hydroperiod was constructed by adding all binary water mask layers and dividing the pixel values by the total number of images, then converting to a relative frequency value between 0 and 100%. The percent frequency at any pixel location, therefore, represents the relative proportion of images that displayed open water at that location and is here assumed analogous to the long-term average seasonal hydroperiod, with 100% representing permanent open water and 0% representing areas of the floodplain that are rarely inundated. Note, 0% is here interpreted as "rarely" inundated, as from the temporal sampling, all that can be determined with confidence is that these areas of the floodplain were never inundated at the times of observation. Consequently, "rarely" inundated areas may be occasionally flooded or never flooded but, in either case, these are the floodplain areas that remain driest for the longest periods.
The hydroperiods for T1 and T2 were then compared on the basis of: (i) relative proportions of time for a given area of open water; and (ii) increases or decreases in permanent open water coverage. To explore the influence of river runoff on hydroperiod during the 35-year study period, regression analysis was performed between the open water area for each individual water mask and the average downstream river discharge observed for the day of image acquisition. As a cursory examination of potential climate change impacts on hydroperiod, river discharge and nearby climate variables of temperature and precipitation were summarized for T1 and T2 and compared to establish if any observed change in hydroperiod was consistent with observed or expected shifts in regional hydro-climatology.
Using the 35-year Landsat TM/OLI hydroperiod as a baseline, the MSS and Radarsat 2 hydroperiods have been compared to evaluate the similarity of information content. Specifically, there are two elements of comparison: (i) the similarity in relative hydroperiod (i.e., temporal frequency) of open water area proportions; (ii) a visual comparison to assess the similarity in floodplain water feature definition.

Long-Term Hydroperiod Comparison
The Landsat TM/OLI floodplain hydroperiods along the stretch of the Upper Columbia River floodplain from Brisco to Spillimacheen for T1 and T2 are presented in Figure 4, with summary statistics in Table 1. Areas of permanent open water (floodplain ponds and river channel) are shown in the darkest shades and labelled with the 100% hydroperiod class. Visibly, the two hydroperiods are similar in terms of discernible major channel and open water pond features but, overall, T2 shows a slightly reduced open water frequency. This is more clearly illustrated in Figure 5, where the comparative hydroperiod open water area (as a percent of total floodplain) has been normalized to 10% hydroperiod quantile intervals in a GIS. (Note: the 0% and 100% hydroperiod represent the absolute observations of "no open water" and "always open water", respectively, while all intermediate observation frequencies have been aggregated to 10% quantile intervals.) From the 20% to 100% observation frequency, T2 demonstrates a consistently reduced area of up to 3% of the floodplain extent ( Figure 5). Meanwhile, T2 shows an almost 8% increase at the 0% observation frequency that represents areas of the floodplain that have not displayed any open water (rarely or never inundated) during the observation record (Table 1)      A reduction in permanent water extent (water bodies present in all observation records) as well as an increase in rarely inundated areas (not inundated within any observation records) from T1 to T2 is discernible when the 100% and 0% hydroperiod areas are isolated and differenced ( Figure 6 and Table 2). Note, by isolating "permanent" and "rarely inundated" in this way, all hydroperiods of 10% to 90% are grouped together as "seasonally inundated". This is necessary, as we have high confidence that the end member regions of 100% and 0% frequencies do, in fact, correspond with continuous A reduction in permanent water extent (water bodies present in all observation records) as well as an increase in rarely inundated areas (not inundated within any observation records) from T1 to T2 is discernible when the 100% and 0% hydroperiod areas are isolated and differenced ( Figure 6 and Table 2). Note, by isolating "permanent" and "rarely inundated" in this way, all hydroperiods of 10% to 90% are grouped together as "seasonally inundated". This is necessary, as we have high confidence that the end member regions of 100% and 0% frequencies do, in fact, correspond with continuous open water or no open water during the image observation samples. However, given the potential for noise in water masks, any pixels showing 10 to 90% have a greater potential to be influenced by omission or commission than the end members. Grouping these intermediate frequencies of open water into a single class therefore mitigates any error that could occur by directly comparing each hydroperiod frequency. Overall, 86.2% of the floodplain showed no change in its permanent, seasonal or non-inundated behaviour, with an 8.4% increase in rarely inundated areas adjacent to the edge of the floodplain and a decrease of 3.5% in permanent open water ( Table 2). (Note, the 3.5% of floodplain area decrease represents a >16% reduction in the original permanent open water extent from T1). Water body shrinkage is most notable for two large ponds on either side of the floodplain and across several smaller ponds (red areas in Figure 6). While very small areas of "rarely inundated" loss (1.4%) and water body expansion (0.5%) are observed, these areas are non-contiguous and most likely represent noise or stochastic changes in floodplain land use and cover.   The deterministic nature of the floodplain hydroperiod is illustrated in Figure 7, where the observed area of each Landsat TM/OLI water mask from 1984 to 2019 (n = 61) is regressed against the average downstream river discharge on the same day. The regression result demonstrates that 87% of the variance in the water mask areas can be explained by Columbia River flow at Nicholson, and this is despite the discharge observation being~65 km downstream of the case study location. Of note, this strong relationship is visibly linear, suggesting that the area of inundation scales linearly with discharge within the range of observations made. Peak flows tend to occur during July (Figure 2), so to ascertain if the relationship altered between the early season rising limb or late season falling limb of the river hydrograph (i.e., identify hysteresis in the relationship), the regression analysis was performed again for all water mask and discharge data up to the end of July (n = 26) and all observations following mid-July (n = 46). (Note, overlap in the sample is necessary to ensure that the regression slope captures both the rising and falling limb of seasonal discharge). In both cases, the regression results continued to be strong (R 2 = 0.93 and 0.85, respectively) and slope remained invariant to within 3 s.f. at 0.035, suggesting that the open water area to discharge relationship holds throughout the hydrologically active season.

Landsat MSS and Radarsat 2 Hydroperiod
The hydroperiod images for Landsat MSS and Radarsat 2 are illustrated in Figure 8 alongside the differences with the long-term Landsat TM/OLI hydroperiod classes for visual comparison. Given the differences in resolution, sample frequency and time periods over which these hydroperiods were generated, they are not expected to be identical. However, a qualitative assessment confirms that all three hydroperiod approaches tend to agree in terms of major large permanent open water bodies, as well as there being some correspondence in many of the rarely or never inundated areas within the floodplain. Two distinct differences are apparent, however: (i) for the SAR hydroperiod image, the area of max or seasonal open water inundation is low and restricted to small buffers around permanent open water bodies, while the area of no observed open water is high (Figure 8, Table 3);  (Figure 7) up to the peak discharge on record results in an area of inundation of 31.56 km 2 . This is 2.90 km 2 greater than the geomorphologically constrained floodplain area, suggesting that at some high discharges, the area to discharge relationship will become non-linear and tend to flatten off, as floodplain water depths and/or flow velocities increase rapidly to compensate for the minimal change in area. The discharge at which the floodplain area of 28.66 km 2 is reached is 662 m 3 s −1 . During the study period, this discharge was only exceeded during the single large flood event in July 2002 (noted above), so it is assumed that this model will hold for most discharges, except those associated with rare low frequency major flood events.

Landsat MSS and Radarsat 2 Hydroperiod
The hydroperiod images for Landsat MSS and Radarsat 2 are illustrated in Figure 8 alongside the differences with the long-term Landsat TM/OLI hydroperiod classes for visual comparison. Given the differences in resolution, sample frequency and time periods over which these hydroperiods were generated, they are not expected to be identical. However, a qualitative assessment confirms that all three hydroperiod approaches tend to agree in terms of major large permanent open water bodies, as well as there being some correspondence in many of the rarely or never inundated areas within the floodplain. Two distinct differences are apparent, however: (i) for the SAR hydroperiod image, the area of max or seasonal open water inundation is low and restricted to small buffers around permanent open water bodies, while the area of no observed open water is high (Figure 8, Table 3); (ii) for MSS, small ponds and river channel water features are not well represented. Notwithstanding the reduced sampling frequency and duration for Radarsat 2 and MSS, the general pattern (or signature) of hydroperiod behaviour across the floodplain is similar for all three methods, with a characteristic "U" shape that signifies that the seasonally inundated area quantile (10 to 90% observation frequencies) individually occupies less floodplain area than either permanent water bodies or rarely inundated areas ( Figure 9). Moreover, the permanent water body areas (Table 3) are close to one another at~4 km 2 and only slightly elevated relative to the temporally comparable TM/OLI permanent water body areas for T1 and T2 (Table 2). However, the systematic differences in open water feature visibility between TM/OLI and both Radarsat 2 and MSS result in an overall reduction in the observed area of inundation across the full hydroperiod frequency range (Figure 9).  (Note, the 0% hydroperiod frequency is not illustrated due to high variability across the three methods; see Tables 2 and 3).

Long-Term Changes in Floodplain Wetland Hydroperiod
Climate change simulations and studies of historical regional hydro-climatology [26][27][28][29][30][31][32][34][35][36]70] suggest tendencies towards regional warming, diminishing snow cover period, increasing rates of snowmelt and evaporation, earlier spring melt peak and lower late summer discharge. Empirical evidence from nearby long-term continuous meteorological (Cranbrook: ~150 km south) and downstream river discharge (Nicholson: ~65 km north-west) records are consistent with these findings and expectations. From Table 4, minimum, mean and maximum temperatures for the more recent T2 period have increased relative to the preceding period, while mean annual total precipitation has decreased. Little can be made of the drop in maximum discharge, which represents two single, very high flow events, but it is noteworthy that the 75th percentile in stream flow has also decreased slightly from T1 to T2, while the overall mean and lower quartile has increased. This observation of increased mean flows is consistent with more rapid snowmelt and potentially elevated glacial melt contributions [26][27][28][29]. However, slightly increased overall river flow volumes may initially appear contradictory to the reduction in hydroperiod ( Figure 5) and the loss of permanent water bodies ( Figure 6) observed in the Landsat TM/OLI record.

Long-Term Changes in Floodplain Wetland Hydroperiod
Climate change simulations and studies of historical regional hydro-climatology [26][27][28][29][30][31][32][34][35][36]70] suggest tendencies towards regional warming, diminishing snow cover period, increasing rates of snowmelt and evaporation, earlier spring melt peak and lower late summer discharge. Empirical evidence from nearby long-term continuous meteorological (Cranbrook:~150 km south) and downstream river discharge (Nicholson:~65 km north-west) records are consistent with these findings and expectations. From Table 4, minimum, mean and maximum temperatures for the more recent T2 period have increased relative to the preceding period, while mean annual total precipitation has decreased. Little can be made of the drop in maximum discharge, which represents two single, very high flow events, but it is noteworthy that the 75th percentile in stream flow has also decreased slightly from T1 to T2, while the overall mean and lower quartile has increased. This observation of increased mean flows is consistent with more rapid snowmelt and potentially elevated glacial melt contributions [26][27][28][29]. However, slightly increased overall river flow volumes may initially appear contradictory to the reduction in hydroperiod ( Figure 5) and the loss of permanent water bodies ( Figure 6) observed in the Landsat TM/OLI record. Table 4. Summaries of daily April to October temperature and total annual precipitation for Cranbrook, and daily Columbia River discharge at Nicholson for T1 (n = 4066 days) and T2 (n = 3638 days). * P25 and P75 refer to the 25th and 75th percentiles (quartiles) of the long-term discharge record. A systematic shift towards higher flows during the spring melt period followed by lower flows in late July to early September is evident in Figure 10. Therefore, while mean flow may have increased and peak flow decreased, it appears that another important factor in the reduced active season hydroperiod behaviour during T2 is the small reduction in river flow during late summer, when evaporation rates are at their highest. A reduction in high summertime flows during late June to early August would inevitably lead to fewer late season overbank events, as well as reduced water exchange between the main river channel and floodplain wetlands via levee breaches and small feeder channels. This reduction in connectivity to runoff in the main channel during late summer means that the isolated and permanent water bodies in the floodplain likely would have been isolated for longer during T2 than T1. Due to the low hydraulic conductivity of floodplain substrates, most loss of water will be due to evapotranspiration. Given that temperatures increased and precipitation decreased in T2 (Table 4), it is reasonable to assume that the potential for floodplain evaporative loss was also greater during T2 than T1. These changes to the seasonal floodplain inundation regime, combined with enhanced losses to late season evaporation, may explain the observed decrease in permananent water bodies from T1 to T2 ( Figure 6, Table 2).
A systematic shift towards higher flows during the spring melt period followed by lower flows in late July to early September is evident in Figure 10. Therefore, while mean flow may have increased and peak flow decreased, it appears that another important factor in the reduced active season hydroperiod behaviour during T2 is the small reduction in river flow during late summer, when evaporation rates are at their highest. A reduction in high summertime flows during late June to early August would inevitably lead to fewer late season overbank events, as well as reduced water exchange between the main river channel and floodplain wetlands via levee breaches and small feeder channels. This reduction in connectivity to runoff in the main channel during late summer means that the isolated and permanent water bodies in the floodplain likely would have been isolated for longer during T2 than T1. Due to the low hydraulic conductivity of floodplain substrates, most loss of water will be due to evapotranspiration. Given that temperatures increased and precipitation decreased in T2 (Table 4), it is reasonable to assume that the potential for floodplain evaporative loss was also greater during T2 than T1. These changes to the seasonal floodplain inundation regime, combined with enhanced losses to late season evaporation, may explain the observed decrease in permananent water bodies from T1 to T2 ( Figure 6, Table 2).

Extending the Hydroperiod
To examine influences on floodplain hydroperiod over a longer time-scale where hydroclimatic changes may be greater or more readily discerned requires the use of lower resolution imagery that pre-dates the Landsat Thematic Mapper mission. The Landsat Multispectral Scanner is an obvious choice due to instrument and acquisition similarities and the potential to extend the record by another 12 years to 1972. However, the reduced spatial resolution of the MSS image pixels (60 m width or

Extending the Hydroperiod
To examine influences on floodplain hydroperiod over a longer time-scale where hydroclimatic changes may be greater or more readily discerned requires the use of lower resolution imagery that pre-dates the Landsat Thematic Mapper mission. The Landsat Multispectral Scanner is an obvious choice due to instrument and acquisition similarities and the potential to extend the record by another 12 years to 1972. However, the reduced spatial resolution of the MSS image pixels (60 m width or 3600 m 2 area) relative to TM/OLI (30 m or 900 m 2 ) means that key pond and river channel water features in the floodplain were not clearly defined or were omitted altogether ( Figure 8). In this part of the Columbia River floodplain, the main channel varies in width between approximately 30 and 40 m, which is why the majority of the main channel is visible as permanent open water in the TM/OLI hydroperiod (Figure 4), yet it either shows up as seasonal or not at all in the MSS hydroperiod ( Figure 8). Consequently, MSS is well suited to hydroperiod generation in areas of slightly larger, more continuous water features, but if employed to extend the Landsat TM/OLI hydroperiod for floodplain areas similar to the Upper Columbia, then some cross-calibration of the hydroperiods over a common time-frame would be recommended.
Continued open water monitoring and hydroperiod construction into the future can be achieved with OLI and many other optical satellite sensors but their primary drawback is cloud cover and the resultant inability to guarantee data on an as-needed basis. With new high spatial and temporal resolution SAR solutions like the Radarsat Constellation Mission [69], reliable repeat open water mapping and monitoring is one of the high-priority mission goals, with automated single polarization (HH) water mask and hydroperiod generation becoming routine tasks [50,51,62,71]. However, it is clear from Figures 8 and 9 and Table 3 that the HH SAR hydroperiod properties are distinct from those from optical sensors, with the area of max or seasonal open water extent being comparatively low and restricted to small buffers around permanent open water bodies, while the area of apparent "dry land" is high. This can be attributed to side-looking radar signals incident upon seasonally inundated water bodies where areas of emergent marshland vegetation will elevate the SAR backscatter intensity above the dB threshold for pure open water. Consequently, Radarsat-derived water masks are unlikely to produce equivalent hydroperiods to those from OLI, due to emergent marsh vegetation appearing "bright" within a slant range single polarisation HH SAR scene, yet largely invisible to nadir pointing optical sensors. Flooded vegetation can be identified using other SAR techniques such as polarimetry, which can help to solve this problem [59,71,72].

Conclusions
Applying a NDVI-MNDWI index subtraction method [55] to a time series of Landsat TM and OLI satellite image data enabled the construction of a continuous hydroperiod map from 1984 to 2019 describing the seasonal and permanent open water characteristics of a section of the Upper Columbia River floodplain. The pattern of inundation over the floodplain wetland complex between Brisco and Spillimacheen was found to be deterministically related to the downstream river flow (R 2 = 0.87, n = 61), a relationship that has clear potential for wetland vulnerability assessment should future runoff simulations become available. The floodplain hydroperiod characteristics of distinct permanent water bodies and rarely inundated areas either side of large areas of seasonal inundation are generally consistent across the TM/OLI early and late time periods, as well as the MSS and Radarsat 2 hydroperiod samples. The most notable divergence between the early (1984 to 2002) and the latter (2003 to 2019) hydroperiod character is the~16% loss of permanent open water (3.5% change in total floodplain area), which corresponds with earlier annual spring melt and peak flows combined with reduced late season flows and enhanced evaporative loss. These observations in the historical record are consistent with expectations for a changing regional climate and suggest that permanent and seasonal open water areas may continue to diminish and thus increase wetland vulnerability. While Radarsat Constellation Mission hydroperiods will be distinct to those derived from optical sources, it is recommended that floodplain wetlands considered vulnerable to changing climates be the subject of high temporal resolution monitoring to track changes in open water hydroperiod and thus provide an early warning of imminent loss of wetland ecosystem services.