The sign, magnitude and potential drivers of change in surface water extent in Canadian tundra

The accelerated rate of warming in the Arctic has considerable implications for all components of ecosystem functioning in the High Northern Latitudes. Changes to hydrological cycle in the Arctic are particularly complex as the observed and projected warming directly impacts permafrost and leads to variable responses in surface water extent which is currently poorly characterized at the regional scale. In this study we take advantage of the 30 plus years of medium resolution (30 m) Landsat data to quantify the spatial patterns of change in the extent of water bodies in the Arctic tundra in Nunavut, Canada. Our results show a divergent pattern of change—growing surface water extent in the north-west and shrinking in the south-east—which is not a function of the overall distribution of surface water in the region. The observed changes cannot be explained by latitudinal stratification, nor is it explained by available temperature and precipitation records. However, the sign of change appears to be consistent within the boundaries of individual watersheds defined by the Canada National Hydro Network based on the random forest analysis. Using land cover maps as a proxy for ecological function we were able to link shrinking tundra water bodies to substrates with shallow soil layers (i.e. bedrock and barren landscapes) with a moderate correlation (R2 = 0.46, p < 0.001). It has previously been reported that rising temperatures are driving a deepening of the active layer and shrinking water bodies can be associated with coarse textured soils beneath the lakes. Unlike water bodies with soil, or gravel, beneath them the water bodies that are situated on bedrock are likely cut off from ground water. Drying water bodies clustered in areas of bedrock and thin soils points to evaporation as an important driver of surface water decrease in these cases.


Introduction
The North American High Northern Latitudes are currently the subject of intensive research focus in both Canada and the United States. Over the past 30 years air temperature in this region has increased faster than in the rest of the world with an average of 2 • C-4 • C (Miller et al 2010, Serreze andBarry 2011) increase over that time span. The observed warming has been identified as the primary driver of a broad spectrum of environmental changes from enhanced vegetation response (i.e. 'greening') (Goetz et al 2010, Ju and Masek 2016, McManus et al 2012, to the deepening of the active layer (seasonally melted layer above permafrost) and talik formation (unfrozen layer of soil beneath lakes and above the permafrost) in the Arctic (Arp et al 2016, Brown et al 1998, Camill 2005, Jorgenson et al 2010, and to an overall trend toward lower sea ice cover (Stroeve et al 2011). Accounting for more than 20% of the land surface across the Arctic tundra, water bodies are a critical component of tundra ecosystems (Downing 2010, Downing et al 2006. These water bodies are essential to the hydrological cycle, carbon cycle and overall energy balance in the region (Bring et al 2016, Cole et al 2007, Francis et al 2009, Stokstad 2004, Walter et al 2007, White et al 2007. Not surprisingly, numerous recent studies have focused on extensive changes across various components of the hydrological cycle including longer inland water ice-free season (Chen et al 2013, Smejkalova et al 2016, increased inland water temperatures (O'Reilly et al 2015, Schneider and Hook 2010, Sharma et al 2015, and a change in the overall extent of surface water across High Northern Latitudes , Jepsen et al 2013, Jorgenson et al 2006, Roach et al 2013, Rover et al 2012, Smith et al 2005, Smol and Douglas 2007b. Many factors including changes in temperature, precipitation, evaporation, vegetation cover, and subsurface soil structure affect the surface water extent of inland water bodies (Bring et al 2016, Gibson and Edwards 2002, Yoshikawa and Hinzman 2003. Incoming precipitation is balanced by runoff, surface storage, evaporation, and subsurface drainage. In high northern latitudes (HNL), however, permafrost prevents surface water from percolating into ground water and thus largely controls subsurface drainage (Jorgenson and Grosse 2016, Jorgenson et al 2010, Kokelj and Jorgenson 2013. In regions of permafrost this, in general, results in water ponding in depressions where there is not enough slope or an outlet to permit overland flow of water. Surface water change in discontinuous permafrost has been shown to be primarily drying due to connection between surface water and sub-surface water. Whereas surface water change in continuous permafrost has been more complex to characterize with the hypothesis of initial wetting followed by drying due to talik formation under the lakes (Smith et al 2005).
Many studies involving the use of satellite imagery have been conducted to ascertain the direction and magnitude of surface water extent change in the HNL over the past 30 years. Most studies of surface water extent change are performed using a thresholding method to classify an image pair (or a few dates) to get the difference in water extent between those dates (Hinkel et al 2007, Smith et al 2005, Yoshikawa and Hinzman 2003. However, those methods that pick a difference between two arbitrarily selected points in time are likely plagued with considerable uncertainty due to a naturally large inter-annual variability of water extent as demonstrated in Carroll and Loboda (2017). Better approaches focus on multi-temporal assessment of long-term trends in surface water change. A time series of Landsat data (one image per year) was used to explore long term change in Yukon Flats, Alaska (Rover et al 2012) which was subsequently used to investigate natural variability in extent of lakes (Chen et al 2013). The most robust currently available methods utilize the full archive of Landsat imagery to account for seasonal variability of surface water extent in addition to the inter-annual variability and obtain the long-term trends in nominal surface water extent as compared to potentially abnormally seasonally high or low gained from a single image per year Loboda 2017, Pekel et al 2016).
The approaches to analyzing the resultant change in surface water extent vary. The global analysis from Pekel et al (2016) compares the difference in a 15 year mean extent of water bodies before and after the year 2000: annual maps from 1984-1999 and 2000-2015 were combined to produce two period maps, respectively, which were then differenced to produce change. Another example represents a regional study with four sites (two in Alaska and two in Siberia), that used Landsat time series with a machine learning approach to perform a per pixel classification into four classes (stable water, stable land, water to land, land to water) based on the trend in spectral indices over time (Nitze et al 2017). Lastly a regional study in northern Nunavut (Canada) used the full time series of Landsat to generate annual maps of the nominal surface water extent. These maps were converted to polygons and used to generate a trend in areal extent for each water body (Carroll and Loboda 2017). Unlike the methods from Carroll and Loboda 2017 and Pekel et al 2016, the method from Nitze et al 2017 produces a single output map that shows the classification with change as a class rather than annual maps that are used in a secondary analysis for change.
Although mapping and quantifying surface water extent change over time at the regional scale is sufficiently challenging, developing an understanding of the drivers of the observed changes at this scale is particularly difficult since datasets representing subsurface parameters are largely unavailable or are grossly oversimplified. The studies that have examined the mechanisms driving surface water extent are limited to local areas with local maps or field measurements (Rover et al 2012) or local airborne data (Jepsen et al 2013, Minsley et al 2012. These studies both concluded that subsurface composition, specifically coarse gravel, was a better determinant of lake change than the overall depth of talik. In a study on Ellesmere Island it was found that lakes on exposed bedrock disappeared completely in the summer of 2006 (Smol and Douglas 2007a). Maps of permafrost (Brown et al 1998) and soils (SLCWG 2010, Tarnocai et al 2002 are coarse in spatial resolution, hence provide little ability to discriminate drivers of change in specific water bodies. Considering the remoteness and difficulty of access for most of circumpolar tundra, large scale studies based on observed (rather than modeled) land surface conditions rely on satellite-derived data products. One approach to developing an understanding of the regional scale drivers is to use other directly observable land surface parameters as a proxy for more general environmental conditions and subsurface characteristics. Specifically, land cover and vegetation fraction-both routinely mapped from satellite observations-can to some degree reflect the type of soils and the depth of active layer in otherwise uniform climatic conditions. Datasets describing land cover (Didiuk and Ferguson 2005, Friedl et al 2010, Olthof and Fraser 2014 and vegetation , Hansen et al 2003, Sexton et al 2013, Townshend et al 2012 are available both globally and regionally. In this paper we use an existing surface water change map (Carroll and Loboda 2017) to explore potential Figure 1. Study region in northern Canada, Nunavut territory where multi-decadal trends in surface water extent have been established (Carroll and Loboda 2017). Water bodies in the region show opposite trends of surface water extent. Red indicates a water body that is increasing in extent while green indicates a water body that is decreasing in surface water extent. There are four weather stations within or just outside of the study region indicated by purple circles. drivers of that change at the regional scale in Nunavut, Canada. Maps of land cover at coarse , Friedl et al 2010, Walker et al 2005 and moderate (Didiuk and Ferguson 2005, Wulder et al 2008 spatial resolution are used to derive proxies for ecological factors from land cover type. The specific objectives of this study are to 1) quantify the spatial relationships and sign (increasing or decreasing) of change shown in Carroll and Loboda 2017, and 2) to relate that change to ecological factors observable at the regional scale including weather station data and various land cover types.

Study area
The study area is located in Nunavut territory primarily in north central continental Canada and extends into the Queen Maud Gulf including parts of several islands (figure 1) bounded by the geographic coordinates 69.2 N, 108.6 W, and 63.5 N, 93.6 W. This area falls within the Arctic Tundra biome (vegetation limited to mosses, grasses, sedges and shrubs) underlain by continuous permafrost (Brown et al 1998). The study area is characterized by low topographic relief and the best available DEM (Canadian Digital Elevation Model used in this study) is at 1:50 000 spatial resolution (∼30 m). The average annual temperature and precipitation range from −11 • C to −14 • C and 273 mm-148 mm from south to north respectively (CCN 2017).
There is limited infrastructure with no major cities or roads in the study domain and the most prominent feature is the Queen Maud Gulf Bird Sanctuary on the northern coast. This sanctuary is the largest federally owned protected area in Canada covering over 61 000 km 2 and is listed in the Ramsar Convention on Wetlands of International Importance. The study area boundaries were selected such that the sanctuary covers approximately half of the continental land area in the study region. This is significant because the protected area designation puts strict limits on the amount of human activity that can be in the protected area, so any changes that are seen here are less likely to be caused by direct human intervention such as dam building or other infrastructure construction.

Data and methods
In our previous work in this study area of North American tundra (Carroll and Loboda 2017) we mapped over 675 000 water bodies, which constitute over 20% of the land cover (excluding ocean), with over 25% of water bodies showing a significant increasing or decreasing trend in surface extent over time (p <0.05). The goal of the current work is to put the previously identified trend in surface water extent into context with landscape and ecology of the region as well as other previously reported changes including the recent findings regarding vegetation greening and browning trends (Goetz et al 2010, Ju and Masek 2016, McManus et al 2012. We conducted a review of available regional-scale land cover and land cover change data to identify key data sets that we used in this analysis (see supplemental section for dataset review).
From our previous work identifying trend in surface water area per water body, we identified the actual areas of observed change (individual per pixel observations) within the water bodies that showed significant increasing or decreasing trends. These grid cells were dichotomized into decreasing and increasing classes by assigning values of −1 for water bodies showing negative trend and 1 for water bodies showing positive trend irrespective of the actual magnitude of the trend. Ocean (Queen Maud Gulf) was masked from all analysis using its maximum extent between 1984-2015 (Carroll and Loboda 2017) to eliminate the coastal shoreline change which is governed by different ecological processes and does not represent the focus of this study. Through the process described above, we identified over 4.3 million grid cells representing ∼3870 km 2 of area with statistically significant trend in long-term surface water extent change within our study region.

Identifying spatial patterns of change
With surface water change identified and quantified down to the grid cell level the next challenge is to identify the driving ecological factors. First order analysis involved identifying spatial patterns of observed dichotomized change. A straightforward way to assess spatial grouping of data is to aggregate from finer resolution to coarser resolution to reveal spatial clustering of areas of change (Lam and Quattrochi 1992). In this case, we took the positive and negative change at 30 m resolution and averaged to 3 km resolution and recorded the percent of grid cells showing change within the resulting 3 km grid cells. Because each 30 m grid cell has a sign (positive or negative) when the averaging occurs, the end result is a 'net' change within the grid cell where equal number of positive and negative values will result in a '0' value within the 3 km grid cell. Spatial aggregation to coarser resolution provides a systematic way to assess the spatial association of nearby change irrespective of other local-scale underlying ecological processes.
Another way to identify spatial patterns is to group the change in terms of other meaningful spatial subsets. We used watershed boundaries to provide a meaningful ecological stratification of the landscape for regionalscale analysis. Surface water area and surface water change were aggregated per watershed, defined by the Canadian National Hydro Network (Natural Resources Canada 2007), using zonal statistics. This produces percent surface water area and percent surface water change per watershed. As in the 3 km percent change, the change per watershed is net change because the sign of change is used in the calculation of area of change.

Random forest analysis
With several datasets available describing surface features it was necessary to determine which datasets are statistically related to surface water change. Random forest analysis (RFA) (Breiman 2001) is a machine learning technique that facilitates assessing the combined influences of multiple inputs on the dependent variable. Unlike many other multivariate statistical models, RFA is not sensitive to the number of independent variables and readily allows for inclusion of both discrete and continuous variables in the analysis. A byproduct of the RFA is a ranking of the independent variables in the context of importance to the RFA. We used this feature of RFA to identify key variables for determining where surface water change will occur. The 4.3 million individual pixels showing change were used as the dependent variable in the RFA. Values for Terrestrial Ecoregions of the World (TEOW), Circum Arctic Vegetation Map (CAVM), Soil Landscapes of Canada (SLCWG), National Hydrography Network-watershed (watersheds), National Hydrography Network-connectivity (connectivity), Canadian DEM-elevation (elevation), Canadian DEM-slope (slope), and MODIS percent bare (bare) were extracted under the raster pixels of change into a table and used as the independent variables in the RFA. A random sample of 10% of the pixels (∼430 000) from the table was used in successive runs of the random forest with increasing numbers of trees starting with five and increasing to 1000 in 30 runs.

Relating spatial patterns of change to ecological drivers
Weather-specifically temperature and precipitationis a critical factor controlling surface water extent with particularly complex conditions in cold regions where ice cover plays a role. Temperature and precipitation data from four weather stations, figure 1, were examined in this study (http://climate.weather.gc.ca/). To establish the overall trend per weather station the daily temperature data was averaged to a single annual average value per year consistent with previous studies (Easterling et al 1997, Jones et al 1999. Most years had between 0 and 7 missing days of data, any year that was missing 10% or more of the observations was set to no data and was excluded from further analysis. Of the four weather stations used only two (Cambridge Bay and Baker Lake) had a consistent record of precipitation while the other two had significant gaps in coverage and were excluded from analysis. The precipitation data (total daily precipitation including frozen precipitation) was summed for each year to produce total annual precipitation for each year from 1985-2015.
There are a number of datasets that describe the land surface of the study region. Many of these land cover characteristics can be used as proxies for ecological parameters. We selected datasets that describe cover type (land cover maps), vegetation condition (Normalized Difference Vegetation Index (NDVI)), elevation and slope, and surface water extent/change across the full study area. A full description of datasets that were investigated is included in the supplemental section. All variables that were used in this analysis provide a continuous representation (i.e. percent, trend, or fraction) except the three classes from the land cover of Canada that describes areas that are either exposed or lichen covered bedrock (henceforth 'barren'). Grid cells that were identified as barren were dichotomized into a separate dataset with two classes-barren and vegetated.
All datasets were aggregated at the watershed level using zonal statistics retaining the average (for continuous) or percent (for discrete) variable per watershed. Univariate statistical relationships between individual surface characteristics and surface water change as the dependent variable were subsequently tested using R 2 values for each variable to determine the goodness of fit.

Results
Initial evaluation of the map of surface water change at 30 m resolution (figure 1) suggests that there are groups or areas within the map where water bodies with positive trend or negative trend are clustered. The clustering becomes more apparent with aggregation of the original map to 3 km grid cells shown in figure 2. The clusters of significant surface water change in the right panel of figure 2 do not simply mimic locations where there are water bodies in the left panel, but rather they show a distinct and different spatial pattern. Visually there does not appear to be a relationship between the amount of surface water area per watershed and surface water change per watershed. The quantitative assessment also shows a poor correlation (R 2 = 0.01) between surface water abundance in the watershed and change. Overall figures 2 and 3 highlight a pattern of change in the study region that goes from north-west to south-east.
Weather is a fundamental driver of surface water change where year to year differences in amount of precipitation can substantially affect the surface water extent for climatologically short periods of time. The purpose of the current study is to identify Figure 3. The study region divided into watersheds using the Canada National Hydro Network. Figure 3(a) on left shows the percent surface water in each of the watersheds. Figure 3(b) on right shows the percent net surface water change per watershed. Both results determined by using the annual surface water extent product (Carroll and Loboda 2017). long term trends in surface water extent and their drivers so we approached analysis of weather data with a similar long term approach. We evaluated temperature and precipitation records from the four official weather stations that are within, or just outside of, the study region with records of at least 30 years (figure 1). Taloyoak, Shepherd Bay and Cambridge Bay stations are located more than 500 km (∼6 • of latitude) to the north of Baker Lake station and thus show substantially lower mean annual temperatures. The annual temperature trend at three of the four stations is positive over the available data record with an average increase of ∼2 • F across the study region (figure 4). This rate of increase is consistent with reports of increasing temperatures throughout the Arctic (Miller et al 2010, Serreze andBarry 2011). The station at Shepherd Bay shows a slight decline in temperatures over the study period. The negative trend at the Shepherd Bay station is strongly influenced by the comparatively high mean annual temperature in the late 1980s. However, a trend beginning from the early 1990s would show a similar direction and magnitude to the trend recorded at all other stations in the study area. Averaged over the four stations, the regional trend in mean annual temperature remains positive ( figure 4(e)).
The annual precipitation, available from only two stations, shows divergent trends for Cambridge Bay (slightly increasing) and Baker Lake (slightly decreasing) (figure 5). Neither trend was significant; however it should be noted that the inter-annual variability was very high at both stations which may indicate an inconsistency in the record, which is not uncommon in cold region precipitation records due to under-catch in winter precipitation (Adam and Lettenmaier 2003). The result of the multivariate analysis of surface parameters as determinants of surface water change from all random forest runs (table 1) shows that the watersheds dataset was the highest ranked predictor of an increasing or decreasing trend in surface water extent followed by elevation, soil types (SLC) and arctic vegetation (CAVM) with an Out of Bag (OOB) error (Breiman 2001) of 48.67%. The bare fraction, terrestrial ecoregions (TEOW), connectivity and slope were the least important in all runs.
Further zonal statistics analysis, performed on the top four variables, was aimed at understanding how these variables are related to the observed pattern of surface water change. The relationship between the percent water change and watersheds can be seen in figure 3(b) with primarily increasing water bodies in the northwest and primarily decreasing water bodies in the southeast. The elevation data provides a per pixel height above sea level that the random forest analysis identified as important on a per pixel level. However, elevation is a continuous variable that, by itself, does not provide a convenient way to identify groups or  spatial patterns without creating arbitrary thresholds of elevation. The soil types (SLC) show only a few broad regions that do not appear to be related to the distinct spatial pattern of surface water change (figure 6). The final predictor from the top four of the random forest analysis is the arctic vegetation map (CAVM) where a single vegetation class-'Cryptogram barren complex (bedrock)'-is linked to a decreasing trend in surface water (figure 7). The spatial resolution (1 km) of CAVM is coarse compared to the features of the region so finer resolution information is needed to confirm and to fully understand the relationship between direction of change and barren or bedrock surface condition. Small outcroppings of bedrock may not be observable in the 1 km base data used to create CAVM and hence the bedrock class may underestimate the true coverage of bedrock in the region.
The Canadian National Hydro Network Watersheds data (Natural Resources Canada 2007), which identify 38 distinct hydrologic units within our study area, are the most ecologically meaningful and the most statistically important (as specified by the random forest analysis) spatial zoning parameter for assessing the relationships between trends in surface water and other satellite-derived surface parameters related to vegetative cover and barren surfaces. Table 2 shows the results of the linear regression between the median of each satellite-derived descriptor and the percent surface water change per watershed. Median value is used to avoid anomalous values (caused by clouds, shadows, or other poor-quality data) that can contaminate an average or maximum calculation.
The percent surface water area per watershed has a low R 2 and p > 0.05 which implies that the amount  and direction of change is not correlated with how much surface water is present in a given watershed. Though the relationship is only moderately strong, the best predictor is the percent bare defined by the barren class of the Landsat ETM+ Land Cover of Northern Canada  with an R 2 of 0.46, figure 8. Though the p value is significant for the analysis of the evapotranspiration (ET) data, the R 2 of 0.15 showing that ET was a weak predictor of surface water change. This is likely due to the coarse spatial resolution of the product, 500 m, and it was only available for the second half of the study period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) because the Terra satellite was not launched until December 1999 and data from the MODIS instrument was not available until February 2000. As a consequence of the spatial resolution it is likely that mixed land and water pixels confound the results. Similarly the Landsat NDVI Vegetation Fraction has a significant p value but a poor R 2 of 0.11 showing it is a weak predictor. In this case it is likely that the NDVI is confounded by soil color or wetness leading to poor correlation.

Discussion
The areal extent of water in a water body at any given time is determined by the hydrological cycle and specifically three primary factors: incoming water (precipitation, inflow from streams/rivers), outgoing water (runoff, evaporation), and topographical depressions that are capable of holding the water. With remotely sensed data we can measure surface features influencing the amount of water stored on that surface. Measurements of subsurface conditions are significantly harder to get because they often require field data collection which can be difficult, expensive, and time consuming. Numerous datasets that describe surface features in various ways are publicly available, however only coarse spatial resolution datasets of subsurface conditions, namely the Circum Arctic Permafrost map (Brown et al 1998) and the Soil Landscapes of Canada (SLCWG 2010), are available at continental scales.
Building upon our previous work (Carroll and Loboda 2017) we have demonstrated that though there is no distinct pattern in fraction of surface water 'cover' there is a pattern from Northwest to Southeast for surface water 'change' in our study region. Direct observations of weather are limited in the study region with only four official weather stations in or near the area that have more than sporadic measurements. Mean annual temperature across the entire region has been increasing at an average rate of ∼0.78 • C per decade over the past 32 years with no distinct spatial variability that could explain the observed spatial patterns of surface water change. While it cannot be ignored that the precipitation data from the two weather stations with complete precipitation data show that precipitation is increasing in the north west (Cambridge Bay) and decreasing in the south east (Baker Lake)-the same pattern that is shown in the surface water changehaving only two stations does not provide sufficient evidence to draw definitive conclusions because there could be significant variation locally.
Watersheds are defined by topographical features that cause water to flow to a specific outlet. The spatial relationship of change determined by location within a watershed shows that there are localized controls on the flow of water. The analysis of surface water change presented here shows that a larger fraction of barren surface and bedrock are associated with the watersheds that decrease in surface water extent. This is significant because the presence of bedrock at or near the surface suggests that the soil is shallow in these locations. Evaluation of surface water change with NDVI (surrogate for vegetation fraction), NDVI trend (change in NDVI over time), and trend in evapotranspiration all yielded weak correlations with surface water change. One possible explanation for this is that further distinction of vegetation type is needed to provide an adequate relationship between surface water change and vegetation (i.e. NDVI is a general measure of vegetation vigor without giving an indication of vegetation type). Generation of a new land cover map to investigate this potential is a large undertaking and well beyond the scope of the work described here.
The observed changes in surface water extent are found throughout the study region. Both air and water temperature Hook 2010, Sharma et al 2015) are increasing in this region. With similar weather and climate conditions for all water bodies (i.e. they do not vary across the study region) the relationship that stands out is the amount of barren or bedrock that is present. In areas where barren or bedrock is less prevalent, hence the vegetation component is higher, the water bodies are expanding. This supports the hypothesis that as the near surface permafrost thaws the water bodies are expanding due to a deepening of the talik beneath the water body. The deeper talik can create a connection to ground water which can have a positive (expansion of water body) or negative (contraction of water body) effect depending on the hydraulic gradient. Talik formation is one possible explanation for the surface water changes that are reported here.
Microtopography is another potential control on the flow of water both on the surface and below through its control on the hydraulic gradient. A DEM finer than the CDED used in this study would be necessary to investigate this and may be possible with future releases of the ArcticDEM (PGC 2017) but at the time of this writing the ArcticDEM is incomplete in the study region which prevents a detailed analysis of the impact of microtopography on surface water. Depth of thaw as well as the composition of the bedrock (fractured or permeable vs solid or impermeable) provide another control on the flow of water which impacts wetting and drying.
The water bodies are drying where barren or bedrock is more prevalent and the vegetation component is low. We hypothesize that the presence of near surface bedrock implies shallow soils in the surrounding area that prevents subsurface lateral water flow and, barring the presence of cracks and fissures in the underlying bedrock, connection to the ground water. There is evidence of similar occurrence in lakes situated over bedrock on Ellesmere Island, Nunavut (Smol and Douglas 2007a) where changes in evaporation to precipitation ratios have resulted in water deficit and hence reduction (or complete loss) of water body extent. In contrast, in areas underlain by deeper soil layers, permafrost thawing associated with rising temperatures creates new pathways for water flow due to melting of ice wedges (Liljedahl et al 2016), deepening talik (Arp et al 2016) or increased surface connectivity between water bodies (Liljedahl et al 2016, Woo andGuan 2006). Change in surface area of water bodies is not uniform and clearly depends on a suite of environmental conditions including subsurface composition, depth of water, as well as evaporation and precipitation ratios.

Conclusions
Over 25% of the 675 000 water bodies in the study region experienced a significant multi-decadal trend in change. Though the entire region experienced the same amount of warming the surface water extent responded differently in different parts of the study domain. The results provided here offer evidence of a direct relationship between surface water change and barren or bedrock substrate at the regional scale. This relationship was identified using two independent datasets (CAVM and Landsat ETM+ Land Cover of Northern Canada) and at different spatial resolutions (1 km and 30 m) respectively. Relating surface water extent change to bedrock detected on the surface provides an extensible method that does not rely on local maps of soils/subsurface or airborne observations with ground penetrating radar. Improved land cover datasets with specific focus on barren (including lichen covered bedrock) land cover types would facilitate expansion of this method to larger region or continental scale. Further improvement would come from a more complete understanding of the subsurface makeup including the amount of fracturing in the bedrock which can also facilitate subsurface flow. providing access to geospatial datasets used in the analysis. We would also like to acknowledge the NASA Center for Climate Studies whose computing support expedited the processing of all of the data. Finally we would like to thank the anonymous reviewers who took the time to review this manuscript.