Multi-Decadal Surface Water Dynamics in North American Tundra

Over the last several decades, warming in the Arctic has outpaced the already impressive increases in global mean temperatures. The impact of these increases in temperature has been observed in a multitude of ecological changes in North American tundra including changes in vegetative cover, depth of active layer, and surface water extent. The low topographic relief and continuous permafrost create an ideal environment for the formation of small water bodies—a definitive feature of tundra surface. In this study, water bodies in Nunavut territory in northern Canada were mapped using a long-term record of remotely sensed observations at 30 m spatial resolution from the Landsat suite of instruments. The temporal trajectories of water extent between 1985 and 2015 were assessed. Over 675,000 water bodies have been identified over the 31-year study period with over 168,000 showing a significant (p < 0.05) trend in surface area. Approximately 55% of water bodies with a significant trend were increasing in size while the remaining 45% were decreasing in size. The overall net trend for water bodies with a significant trend is 0.009 ha year−1 per water body.


Introduction
The North American tundra is a complex landscape where vegetation is interspersed with over one million water bodies, the majority of which are much smaller than 100 hectares (ha) [1,2].Understanding where the water is located and how it is changing over time is a critical component of understanding the role of water in the carbon cycle [3,4], albedo [5,6], energy balance [7][8][9], and quantifying habitat for migratory wildlife [10,11].In the context of global warming, the Arctic region is warming much faster than lower latitudes [12][13][14].The Arctic also has a greater proportion of terrestrial surface occupied by surface water than lower latitudes [2,15,16].The sediments in Arctic lakes are known to have high concentrations of carbon that can be released to the atmosphere if the surface water changes and/or if the active layer deepens [17][18][19][20][21].Given the higher proportion of surface water in the Arctic, the large amount of carbon in the terrestrial Arctic, and the observations of changes in surface water it is important to understand not only the extent of surface water in the Arctic but also how that surface water is changing through time.
In recent years, the extent of surface water has been mapped with several approaches and at regional to global scales.Each of these products has represented an advance in the mapping capabilities but also still has significant limitations including spatial resolution [22][23][24], temporal resolution [25,26] or both [27].Most recently, work has been done to generate continental scale maps of surface water extent at the decadal time step using a time series of Landsat data as input [28,29].Unlike many other Remote Sens. 2017, 9, 497 2 of 15 30 m spatial resolution water maps, the decadal water maps use three years of inputs to generate a single map in order to have enough repeated observations throughout the study region to have confidence that the water that was detected was not a false detection caused by cloud shadows or burn scars, and did not represent a local minimum or maximum caused by short-term weather events (i.e., drought or flood).These decadal water maps represent an advance in the capability to determine surface water change but still do not provide the temporal detail necessary to determine if these decadal maps represent actual long-term trends in water dynamics.
In addition to the mapping efforts there have been attempts to quantify change in surface water extent over time.Most of these change studies [1,30,31] quantify change between two distinct time periods and focus on either a few distinct lakes or the overall regional change.This methodology works reasonably well for mapping forest extent or urban growth that are distinct and comparatively stable, but is less effective with water because water bodies are highly dynamic both inter-annually and seasonally.Quantifying change in surface water extent over large regional to continental areas requires a different approach.
There are two main challenges with generating accurate maps of surface water change.First is having an accurate map of the nominal (not maximum or minimum) extent of water.Second is having enough representative maps in a time series to ensure that any detected change is an actual long term change and not simply seasonal or inter-annual fluctuations in the size of the water body.This can become more problematic when the water bodies are small and shallow because they are more susceptible to short term weather effects (heavy rains or long dry spells).
Inland water accounts for between 10% and 20% of the continental surface in northern Canada depending on the source of the water map [22,23,26].Previous work focused on the abundance [2,22] and importance of small water bodies [32] both globally and in the Arctic.Though the overall coverage of surface water is significant the impact is often overlooked in regional/global models that only consider large contiguous water bodies [9,33].Even studies of volumetric storage and movement of water do not completely include small water bodies [8,[34][35][36].Here, the focus is on generating a series of maps for the purpose of quantifying change rather than simply to map the location and extent of surface water.
The primary goal of this work is to determine the extent of surface water change over a representative region of the North American tundra in north central Canada.To accomplish this goal the current study aims to: (1) create annual maps of nominal surface water extent (1985-2015); (2) quantify the trajectory of overall change; and (3) assess and quantify the trajectory and trend in surface area for individual water bodies.For this work, water bodies are defined as all open surface water including lakes and rivers but excluding oceans and wetlands dominated by emergent vegetation.

Study Area
This research will focus on an area in northern Nunavut territory of Canada (Figure 1) that is bounded by a region with upper left 105 • W, 70 • N and lower right 95 • W, 65 • N.This area is in the Southern Arctic ecoregion [37] which is characterized by wide expanses of shrublands, wet sedge meadows, low topography and small polygonal lakes.The area is underlain by continuous permafrost [38], soils are relatively consistent throughout the area [39] though the best available dataset is very coarse (1:10,000,000), it is away from most human settlements, and away from most industry (oil, gas, and other mining) [40] all of which could be factors in a land cover analysis.These factors, if present, could influence water levels through human extraction of water from lakes or diversion of water from rivers and streams for consumption or industrial use.
Finally, the study area includes a large protected area, the Queen Maud Gulf Bird sanctuary, which also decreases the likelihood of disruptive human activity affecting the landscape.

Methods
The Landsat data archive over Canada is dense (both temporally and spatially) and is publicly available through the US Geologic Survey [41].The Dynamic Surface Water Extent (DSWE) algorithm [42] uses the Landsat surface reflectance data to produce maps of surface water extent for each Landsat scene that is available in the archive.This results in a standardized mapping of water across the Landsat suite of instruments including Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+).Using the full suite of Landsat data enables the creation of a time series of maps since 1984.Annual maps of water extent can provide a rich time series to evaluate change in surface water extent over the past 31 years of the satellite record.
Mapping water in cold regions using visible, near infra-red, and short wave infra-red (such as the Landsat instruments) limits the inputs to the ice-free season, defined here as June through September.Each of the Landsat satellites images every place on Earth every 16 days globally but more frequently at high northern latitudes due to the orbital overlap that increases near the poles [43].This results in a maximum possible observations of ~35 daytime observations over the study region per year.Given that cloud cover is frequent in the region [44] and the persistence of ice in the centers of lakes, the number of potential cloud free images is reduced to 3-8 cloud/ice free observations per year, based on evaluation of data in the current study.To ensure that there are enough observations for the maps to depict nominal (neither maximum nor minimum) extent, a rolling set of cloud free observations over three years of inputs are used to create a single map (e.g., data from 1984, 1985, 1986 used for 1985 map; data from 1985, 1986, 1987 used for 1986 map; etc.).This approach was used effectively with Landsat data to generate the decadal water maps in the high northern latitudes [28].
For the current study, the DSWE data product [42] was selected as input for the annual water maps instead of the custom classification algorithm in the decadal water maps.The DSWE algorithm The study area is located primarily in the Southern Arctic ecoregion and is characterized by low topography and numerous small to moderate sized water bodies.The Queen Maud Gulf Bird sanctuary is indicated with the red polygon.

Methods
The Landsat data archive over Canada is dense (both temporally and spatially) and is publicly available through the US Geologic Survey [41].The Dynamic Surface Water Extent (DSWE) algorithm [42] uses the Landsat surface reflectance data to produce maps of surface water extent for each Landsat scene that is available in the archive.This results in a standardized mapping of water across the Landsat suite of instruments including Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+).Using the full suite of Landsat data enables the creation of a time series of maps since 1984.Annual maps of water extent can provide a rich time series to evaluate change in surface water extent over the past 31 years of the satellite record.
Mapping water in cold regions using visible, near infra-red, and short wave infra-red (such as the Landsat instruments) limits the inputs to the ice-free season, defined here as June through September.Each of the Landsat satellites images every place on Earth every 16 days globally but more frequently at high northern latitudes due to the orbital overlap that increases near the poles [43].This results in a maximum possible observations of ~35 daytime observations over the study region per year.Given that cloud cover is frequent in the region [44] and the persistence of ice in the centers of lakes, the number of potential cloud free images is reduced to 3-8 cloud/ice free observations per year, based on evaluation of data in the current study.To ensure that there are enough observations for the maps to depict nominal (neither maximum nor minimum) extent, a rolling set of cloud free observations over three years of inputs are used to create a single map (e.g., data from 1984, 1985, 1986 used for 1985 map; data from 1985, 1986, 1987 used for 1986 map; etc.).This approach was used effectively with Landsat data to generate the decadal water maps in the high northern latitudes [28].
For the current study, the DSWE data product [42] was selected as input for the annual water maps instead of the custom classification algorithm in the decadal water maps.The DSWE algorithm uses a hierarchical series of spectral tests on the visible and short wave infrared bands of Landsat to determine the presence of water.This algorithm is run on Landsat surface reflectance data that have been ortho-rectified and terrain corrected.The ortho-rectification has a stated error of approximately 1   2   pixel for over 96% of the pixels [45].The Landsat data (from Landsat 5 and 7) used in this study are provided in standard Landsat WRS-2 grid (path 35-42 row 11-15).In this format, there is a significant amount of overlap between adjacent scenes that increases the total number of observations but is not conducive to large area studies.For this study, the DSWE data were processed into annual maps, projected and mosaicked into Canada Alber's Equal Area projection, full description to follow.

Annual Product Generation
The "Raw" unmasked DSWE product was used as input for the annual map generation.The DSWE product has four output values: A combination of high and low confidence water was processed through the same post classification procedure (described below) that was used to generate the decadal water maps [28] to produce annual maps (Figure 2).Though the partial water class correctly identifies edges of water features it also seems to identify many false positives and hence was treated as "not water" in the annual map creation.
Remote Sens. 2017, 9, 497 4 of 15 uses a hierarchical series of spectral tests on the visible and short wave infrared bands of Landsat to determine the presence of water.This algorithm is run on Landsat surface reflectance data that have been ortho-rectified and terrain corrected.The ortho-rectification has a stated error of approximately ½ pixel for over 96% of the pixels [45].The Landsat data (from Landsat 5 and 7) used in this study are provided in standard Landsat WRS-2 grid (path 35-42 row 11-15).In this format, there is a significant amount of overlap between adjacent scenes that increases the total number of observations but is not conducive to large area studies.For this study, the DSWE data were processed into annual maps, projected and mosaicked into Canada Alber's Equal Area projection, full description to follow.

Annual Product Generation
The "Raw" unmasked DSWE product was used as input for the annual map generation.The DSWE product has four output values: A combination of high and low confidence water was processed through the same post classification procedure (described below) that was used to generate the decadal water maps [28] to produce annual maps (Figure 2).Though the partial water class correctly identifies edges of water features it also seems to identify many false positives and hence was treated as "not water" in the annual map creation.Binary maps of water/not water were made from DSWE for each time period.Although each time period included the total of three years, more than one DSWE layer was available for each year.Moreover, the clear surface views within each scene are impacted by clouds and shadows resulting in a varying number of total valid observations for each pixel within the scene.These were summed to get "total water" and "total land" for the three year input period.Probability of water was calculated per pixel as where Probw is the probability of a pixel being water, Wt is the total observations of water for the period, and Lt is total observations of land for the period.Water in the annual map is anything that has greater than or equal to 50% probability of being water.Following this methodology, annual maps were generated for each year from 1985 to 2015; full description can be found in [28].Binary maps of water/not water were made from DSWE for each time period.Although each time period included the total of three years, more than one DSWE layer was available for each year.Moreover, the clear surface views within each scene are impacted by clouds and shadows resulting in a varying number of total valid observations for each pixel within the scene.These were summed to get "total water" and "total land" for the three year input period.Probability of water was calculated per pixel as Prob w = (W t /(W t + L t )) where Prob w is the probability of a pixel being water, W t is the total observations of water for the period, and L t is total observations of land for the period.Water in the annual map is anything that has greater than or equal to 50% probability of being water.Following this methodology, annual maps were generated for each year from 1985 to 2015; full description can be found in [28].

Annual Product Accuracy Assessment
Accuracy assessment was performed by comparing our delineations with very high resolution (VHR) data following the current best practices for remotely sensed data [46].Commercial VHR data are provided to NASA funded scientists through a contract with the National Geospatial Intelligence Agency [47].Multi-spectral data from WorldView-2 were identified for 10 areas within the study region (Figure 3).WorldView-2 multi-spectral data have a nominal spatial resolution of 2 m and provide data in up to 8 bands (4 used in this study: Red, Green, Blue, Near-Infrared) [48].These data were ortho-rectified using the Ames Stereo Pipeline open source software [49].A random sample of 643 points (see Olofsson et al. [46] for procedures for sample size determination) within the boundaries of the VHR data in the study area was chosen for evaluation.Each point was inspected by visual inspection and assigned a value of water or land.

Annual Product Accuracy Assessment
Accuracy assessment was performed by comparing our delineations with very high resolution (VHR) data following the current best practices for remotely sensed data [46].Commercial VHR data are provided to NASA funded scientists through a contract with the National Geospatial Intelligence Agency [47].Multi-spectral data from WorldView-2 were identified for 10 areas within the study region (Figure 3).WorldView-2 multi-spectral data have a nominal spatial resolution of 2 m and provide data in up to 8 bands (4 used in this study: Red, Green, Blue, Near-Infrared) [48].These data were ortho-rectified using the Ames Stereo Pipeline open source software [49].A random sample of 643 points (see Olofsson et al. [46] for procedures for sample size determination) within the boundaries of the VHR data in the study area was chosen for evaluation.Each point was inspected by visual inspection and assigned a value of water or land.These results were compared to the 2010 annual water map because most of the VHR data used in the analysis were acquired in this year.

Identifying Unique Water Bodies
The overall goal of this analysis is to determine if a trend is present for individual water bodies in the study region, thus requiring object-oriented analysis (i.e., changes in discrete water bodies).We define a "water body" as any group of adjacent pixels and water body is inclusive of lakes, rivers and ocean.Each annual raster map was vectorized and individual water body objects were defined using the "queen" adjacency rule where any pixel that touches the candidate pixel is included in the polygon and the area of the polygon was calculated.However, the highly dynamic nature of small and shallow water bodies determines that the water bodies can split, merge, disappear, or appear over time.To ensure that these changes are captured in our analysis, a master map was generated These results were compared to the 2010 annual water map because most of the VHR data used in the analysis were acquired in this year.

Identifying Unique Water Bodies
The overall goal of this analysis is to determine if a trend is present for individual water bodies in the study region, thus requiring object-oriented analysis (i.e., changes in discrete water bodies).We define a "water body" as any group of adjacent pixels and water body is inclusive of lakes, rivers and ocean.Each annual raster map was vectorized and individual water body objects were defined using the "queen" adjacency rule where any pixel that touches the candidate pixel is included in the polygon and the area of the polygon was calculated.However, the highly dynamic nature of small and shallow water bodies determines that the water bodies can split, merge, disappear, or appear over time.To ensure that these changes are captured in our analysis, a master map was generated showing the extent of water over the entire 31-year study period.
For this project, annual maps were created that show land and water as 0 and 1, respectively.Calculating the sum of these 31 maps yields a single map showing the frequency of occurrence of water per pixel giving values from 0-31.Using the adjacency rules described earlier, we identify all contiguous water bodies with a frequency of occurrence of 1-31.The resulting master map is the maximum possible extent of a water body over the 31-year study period.In an individual year it is possible for there to be more than one distinguishable water body but we want to consider them as part of the same water body through time to make sense of the area statistics (Figure 4).The master map was used to specify water body object identifiers which were propagated through the full temporal extent of the data record.Subsequently, area of individual water body objects was calculated at the annual time step and the total change in water extent is quantified on the per object basis.
Remote Sens. 2017, 9, 497 6 of 15 For this project, annual maps were created that show land and water as 0 and 1, respectively.Calculating the sum of these 31 maps yields a single map showing the frequency of occurrence of water per pixel giving values from 0-31.Using the adjacency rules described earlier, we identify all contiguous water bodies with a frequency of occurrence of 1-31.The resulting master map is the maximum possible extent of a water body over the 31-year study period.In an individual year it is possible for there to be more than one distinguishable water body but we want to consider them as part of the same water body through time to make sense of the area statistics (Figure 4).The master map was used to specify water body object identifiers which were propagated through the full temporal extent of the data record.Subsequently, area of individual water body objects was calculated at the annual time step and the total change in water extent is quantified on the per object basis.The combined record of annual area estimates for each of the water body objects was processed in R statistical analysis software to generate a linear regression (ordinary least squares regression) through time.The slope, intercept, correlation coefficient (r 2 ), Root Mean Square Error (RMSE), and probability-value (p-value) were calculated and recorded.

Accuracy Assessment Results
The results of the accuracy assessment show that the annual water maps for year 2010 compared to WorldView-2 data represent water distribution of the tundra landscape of Northern Canada with the overall accuracy of 95% (Table 1).Although the producer's accuracy for water bodies was lower (87%), visual inspection of the erroneous pixels shows that the misclassified pixels occurred at the The combined record of annual area estimates for each of the water body objects was processed in R statistical analysis software to generate a linear regression (ordinary least squares regression) through time.The slope, intercept, correlation coefficient (r 2 ), Root Mean Square Error (RMSE), and probability-value (p-value) were calculated and recorded.

Accuracy Assessment Results
The results of the accuracy assessment show that the annual water maps for year 2010 compared to WorldView-2 data represent water distribution of the tundra landscape of Northern Canada with the overall accuracy of 95% (Table 1).Although the producer's accuracy for water bodies was lower (87%), visual inspection of the erroneous pixels shows that the misclassified pixels occurred at the edge of a water body or in a pond smaller than one Landsat pixel.

Long Term Water Dynamics
The study region encompasses over 30,000,000 ha divided into two categories: land 25,820,000 ha (84.4%) and inland water 4,784,100 ha (15.6%), all values based on average area of our maps over the 31-year study period.Oceans were eliminated from evaluation as were ocean coastlines because changes to oceans and coastlines are outside the scope of the study question.Over 675,000 individual water bodies were identified during the 31-year study period, using the master map of contiguous water bodies described earlier.The inland water bodies range in size from 0.09 ha (one Landsat pixel) to 350,000 ha based on average size over the 31-year study period (Table 2).The majority of detected water bodies (67%) are small (<1 ha) and the total number of water bodies decreases as size class increases.
Comparing maps from individual years exposes differences in extent of water bodies between those years.The total area of inland water varies annually from a low in 1999 of 4,435,692 ha to a high of 4,764,440 ha in 1989 (Figure 5) with a significant amount of inter-annual variability.The three red circles are the area values at 10-year intervals beginning with 1989 suggesting a local maximum, minimum and another maximum.Water bodies that are less than 1 ha routinely fluctuate in size through the data record.By overlaying multiple years together, it is possible to see the progression of change in specific water bodies (Figure 6).Typical fluvial processes can be seen as rivers meander and midstream islands change shape over time (Figure 7).The amount and direction of change was quantified by performing a linear regression analysis (ordinary least squares) on area of water bodies over 31 years and is reported in Table 3.Nearly 87% of all individual water bodies in the study region showed a trend in area with 45% showing positive trend (increasing size) and 42% showing negative trend (decreasing size).Over 168,000 (25%) water bodies that show a trend had p-value < 0.05 and nearly 73,000 (11%) had p-value < 0.01.Small (0.1-1 ha) water bodies constitute the majority of water objects with significant (at p < 0.05) trend in surface water area (Table 4).This distribution is similar to the overall size distribution of water bodies (Table 2).The results show that for the majority of water bodies, specifically those between 0.1 and 1000 ha in size, the fraction of bodies that exhibited a significant trend in surface water area at p < 0.05 is approximately 27%.This fraction is smaller for very small (<0.1 ha) bodies at 21% and large water bodies at 24%, 14%, and 22% for 1000-10,000 ha, 10,000-100,000 ha, and >100,000 ha, respectively.
Breaking this down further into water bodies that are increasing compared to those that are decreasing more differences in the trends become apparent.The amount and direction of change was quantified by performing a linear regression analysis (ordinary least squares) on area of water bodies over 31 years and is reported in Table 3.Nearly 87% of all individual water bodies in the study region showed a trend in area with 45% showing positive trend (increasing size) and 42% showing negative trend (decreasing size).Over 168,000 (25%) water bodies that show a trend had p-value < 0.05 and nearly 73,000 (11%) had p-value < 0.01.The smallest size class shows a distinct trend toward decreasing in size: 61% of all water bodies in size <0.1 ha that exhibit a significant trend decreased in extent over time.Since a single Landsat pixel is ~0.09 ha any water body that is <0.1 ha that decreased in size either disappeared or became too small to be detected by Landsat anymore.In comparison, the next three size classes show a distinct trend toward increasing in size at approximately the same proportion (61-63% of water bodies that exhibit a significant trend in each category).Approximately half of the medium (52% of the 1000-10,000 ha) and most large (75% and 100% of 10,000-100,000 ha and >100,000 ha, respectively) water bodies with significant trend have decreased in size over time.The analysis here is possible because the use of the master map prevents the counts from inflating artificially when several water bodies coalesce into a single water body or conversely a single water body dries and splits into two or more water bodies.The master map keeps those water bodies together throughout the analysis.The water bodies that show a significant change can be found throughout the study area (Figure 8).Water bodies that exhibit a trend toward growth are colored green, those that show a trend towards shrinking are colored red and those water bodies that have no significant trend are light blue-grey.There are water bodies that are increasing in size that are spatially adjacent to water bodies that are decreasing in size.It is also apparent that the water bodies are distributed throughout the study region regardless of size, i.e., small lakes and large lakes can be found anywhere in the study region.

Discussion
In this study, the annual maps are generated from a rolling set of three years of inputs which results in a map of nominal extent of surface water over the period of the inputs.By doing this the likelihood that an individual map is showing maximum or minimum extent is reduced thereby increasing the validity of inter-comparison with other maps produced in the same way, i.e., the time series created here.Accuracy assessment was performed on the map for the year 2010 to closely

Discussion
In this study, the annual maps are generated from a rolling set of three years of inputs which results in a map of nominal extent of surface water over the period of the inputs.By doing this the likelihood that an individual map is showing maximum or minimum extent is reduced thereby increasing the validity of inter-comparison with other maps produced in the same way, i.e., the time series created here.Accuracy assessment was performed on the map for the year 2010 to closely match the VHR data that were used for the assessment.It is assumed that the accuracy is similar for maps from other years because the input product (DSWE) and the methodology remain the same.
Figure 4 suggests that the overall area of surface water changes from year to year.This clearly shows that studies that use any two years as inputs and employ simple differencing to identify change will inevitably produce erroneous results.The majority of water bodies that show significant trends are growing in size at an average rate of 0.03 ha per year.The net trend in area of surface water per water body is 0.009 ha calculated as average trend in growth minus average trend in decline.This finding is consistent with predictions that as the Arctic warms the active layer of permafrost deepens [50,51].As the active layer deepens, this can result in the erosion of the margins of a water body often resulting in expansion of the water body.Conversely, deepening of the active layer could result in the collapse of the lake bottom which would likely cause the draining of the water body [10].The changes to a water body by deepening of the active layer are not necessarily uni-directional; in some cases, a water body could initially expand until a critical threshold is exceeded at which point an edge fails or the bottom falls out and the water body drains (i.e., thermokarst).
The novelty in this analysis comes in the object based analysis that tracks water body area through time by spatial location.Through this analysis any water body that is shown to disappear has not been absorbed by the growth of a larger immediately adjacent water body.Table 4 shows that size of water body plays a role in the direction of change (increasing or decreasing) where the smallest water bodies (<0.01 ha) are decreasing in size over the study period and the remaining small water bodies (<1 ha) are increasing in size.This could imply that water from the smallest water bodies is moving laterally into larger adjacent water bodies or it is exiting the system in a different way (possibly evaporation or drainage into ground water).However, determining exactly where the water is going is beyond the scope of this study.
While there are specific instances of complete drainage of a water body or the creation of new water bodies, these instances are rare.The more interesting result is the quantification of a significant trend showing expansion of surface water area over the long time series.While the changes are distributed throughout the study region (Figure 8) there appear to be regions of increasing and decreasing water bodies.The water bodies that are increasing appear to be concentrated in the northwest portion of the image and the water bodies that are decreasing appear to be concentrated on the eastern half of the image.Further investigation with ancillary data such as ecoregions, soils, elevation and vegetation cover may provide additional insight into the potential patterns of change.The implication of wide spread growth in surface water lies in the extensive stores of carbon in the soils and sediments in the region.Greater surface water area can lead to (or be caused by) deepening of the active layer beneath and around the water bodies.One potential implication of changing surface area of water is the possibility of mobilizing the significant amount of soil based carbon to the atmosphere through methanogenesis, denitritrification, aerobic decomposition/respiration and other biogeochemical processes.
The methodology used to create the annual water maps from a series of water classifications is robust and has been employed for both Landsat and MODIS instruments [22,28] and is similar to the method adopted for the Global Landsat Water maps [27].The development of these new time series products of surface water make it possible to apply our object based analysis methodology in other regions.The challenge arises in large regions with a rapid increase in the number of objects.The current process with over 600,000 objects pushed the limits of what a standard software package (ArcGIS in this study) would handle.Significant improvements in data handling in this software or the development of custom routines that are outside of proprietary software are needed to facilitate processing larger numbers of objects.Ideally, these new processes would be developed to take advantage of "Big Data" processing environments and parallel processing to both speed up the processing time and enable larger datasets to be considered.
The DSWE water detection algorithm was originally designed and tuned to detect water in the contiguous United States.This is the first application of the algorithm in the tundra ecotone and the algorithm has performed well as demonstrated by the validation presented in this text.Future enhancements to the DSWE algorithm including improvements to the "partial water" class will have important implications for advancing the methodology presented here to support studies that quantify carbon exchange in the tundra.In addition, the identification of a snow/ice coupled with improved cloud detection would increase the overall number of observations that could be used.These improvements have been suggested to the developers of the DSWE and are currently under consideration.
Previous analyses of water body change in the Arctic have focused on only a few observations and there has been more work done in the Alaskan North Slope than the Canadian high Arctic [52][53][54][55].Coarse resolution analysis using MODIS data [1] was limited to differencing over a decade and at 250 m spatial resolution does not fully capture the dynamics because so many of the water bodies are well below detection levels with MODIS (1 MODIS pixel at 250 m resolution ~6.25 hectares) see Table 4 showing over half of the water bodies with a significant trend are 1 hectare or less.Coarser resolution analysis (25 km) using microwave data [56] has the advantages of fractional coverage and the ability to collect data under cloudy conditions but is still too coarse to capture the dynamics of small water bodies.Most recently global annual water maps have been generated from Landsat [27], however the change mapping provided combines all maps up to 1999, all maps after 1999 (yielding two maps) and then does a simple difference.Based on the information in Figure 5, this approach is likely to lead to erroneous results at least for the current study region.The ability to use the full time series of Landsat data to generate a time series of object based (individual water bodies) measurements provides a new perspective on the true dynamics of water bodies in the region.The density and diversity of remote sensing observations has increased with the launch of Sentinel 1A/B (radar observations) and Sentinel 2 A/B (10-20 m visible and near infrared observations).Combining the historical time series with the current and potential future observations from the Sentinel constellation will provide an extended record into the future.

Conclusions
This is one of the first analyses to produce annual water maps that show the nominal water extent in each year over 31 years for a region.A large inter-annual variability in the results demonstrates that it is necessary to do time series analysis to fully understand the dynamics of surface water in the region.Size analysis of water bodies confirms previous studies that show that small water bodies dominate the landscape.Our results show that over the past 31 years the smallest water bodies have been shrinking or disappearing entirely, whereas slightly larger to moderate water bodies have been increasing in size.The area based linear regression analysis shows an overall net growth of 0.009 ha year −1 per water body.The occurrence of significant trend in surface water area in the study area, which is located far from most anthropogenic activity, implies that environmental or climate factors are driving the observed change.Further investigation is required to determine if the potential periodicity seen in the graph of total area of surface water over 31 years (i.e., Figure 5) is related to a climate driver or other environmental factors and an analysis of spatial patterns and clustering will reveal potential areas for further investigation.

Figure 1 .
Figure 1.Study area in North American Arctic region, north central Nunavut territory Canada.The study area is located primarily in the Southern Arctic ecoregion and is characterized by low topography and numerous small to moderate sized water bodies.The Queen Maud Gulf Bird sanctuary is indicated with the red polygon.

Figure 1 .
Figure 1.Study area in North American Arctic region, north central Nunavut territory Canada.The study area is located primarily in the Southern Arctic ecoregion and is characterized by low topography and numerous small to moderate sized water bodies.The Queen Maud Gulf Bird sanctuary is indicated with the red polygon.

Figure 2 .
Figure 2. Algorithm flow for the generation of annual water maps from the individual dates of DSWE (reproduced from [28]).In Panel (a) (left most panel), individual scenes are converted from four classes to two (land and water) then summed to get total observations of land and total observations of water for the period.Panel (b) shows the "total water" for an individual path/row.Panel (c) shows the mosaicked non-overlapping path/rows which are then summed.Panel (d) shows the final "total water" for the full region of interest.

Figure 2 .
Figure 2. Algorithm flow for the generation of annual water maps from the individual dates of DSWE (reproduced from [28]).In Panel (a) (left most panel), individual scenes are converted from four classes to two (land and water) then summed to get total observations of land and total observations of water for the period.Panel (b) shows the "total water" for an individual path/row.Panel (c) shows the mosaicked non-overlapping path/rows which are then summed.Panel (d) shows the final "total water" for the full region of interest.

Figure 3 .
Figure 3. Distribution of WorldView-2 (WV2) scenes used in the accuracy assessment of the annual water maps.Footprints of WV2 scenes are shown as dark grey rectangles distributed throughout the image.

Figure 3 .
Figure 3. Distribution of WorldView-2 (WV2) scenes used in the accuracy assessment of the annual water maps.Footprints of WV2 scenes are shown as dark grey rectangles distributed throughout the image.

Figure 4 .
Figure 4.In the images above water is shown in black and land is shown in light grey.The first five images show a lake complex (small unnamed lake in northern Nunavut) in five individual years.The final image shows the master map which is the maximum extent in the whole 31-year record.Through time you see the lake shrinks and splits into several components.The master map allows all of these components to be related to the same water body even when they split off into individual pieces.The years shown here are chosen as representative examples of the 31-year record.

Figure 4 .
Figure 4.In the images above water is shown in black and land is shown in light grey.The first five images show a lake complex (small unnamed lake in northern Nunavut) in five individual years.The final image shows the master map which is the maximum extent in the whole 31-year record.Through time you see the lake shrinks and splits into several components.The master map allows all of these components to be related to the same water body even when they split off into individual pieces.The years shown here are chosen as representative examples of the 31-year record.

Figure 5 .
Figure 5.Total annual area of surface water between 1985 and 2015.Red circles denote local temporal maxima and the minimum in the record.Figure 5. Total annual area of surface water between 1985 and 2015.Red circles denote local temporal maxima and the minimum in the record.

Figure 5 . 15 Figure 6 .
Figure 5.Total annual area of surface water between 1985 and 2015.Red circles denote local temporal maxima and the minimum in the record.Figure 5. Total annual area of surface water between 1985 and 2015.Red circles denote local temporal maxima and the minimum in the record.Remote Sens. 2017, 9, 497 8 of 15

Figure 6 .
Figure 6.Difference in extent for several lakes in the study region.The lighter colors indicate that water was present in early years of the study but not present in later years.

Figure 6 .
Figure 6.Difference in extent for several lakes in the study region.The lighter colors indicate that water was present in early years of the study but not present in later years.

Figure 7 .
Figure 7. Variability of river extent and flow from 1985-2015.Lighter colors indicate that water was present in some years but not in all years.This is particularly noticeable in the edges of the river and in the islands in the middle of the channel.

Figure 7 .
Figure 7. Variability of river extent and flow from 1985-2015.Lighter colors indicate that water was present in some years but not in all years.This is particularly noticeable in the edges of the river and in the islands in the middle of the channel.

Figure 8 .
Figure 8. Spatial distribution of water bodies with a significant (p < 0.05) trend in surface water area (ha/year) over the 31-year time period of the study.Water bodies that are increasing in size are shown in green while the ones that are decreasing in size are shown in red.

Figure 8 .
Figure 8. Spatial distribution of water bodies with a significant (p < 0.05) trend in surface water area (ha/year) over the 31-year time period of the study.Water bodies that are increasing in size are shown in green while the ones that are decreasing in size are shown in red.

Table 1 .
Confusion matrix for accuracy assessment of annual water maps using WorldView-2 multi-spectral data.

Table 2 .
Distribution of water bodies detected in the study region by size.All sizes are in hectares and based on average area over the 31-year study period (the very low fractional representation of the largest water bodies in the total water body count necessitates the use of precision values to the third decimal point).
largest water bodies in the total water body count necessitates the use of precision values to the third decimal point).

Table 3 .
Overall 31-year trend in size of water bodies.

Table 4 .
Distribution of water bodies that exhibit significant trend in surface water area at p < 0.05 by size and trajectory of change.