The Canadian Water Resource Vulnerability Index to Permafrost Thaw (CWRVI_PT)

Vulnerability is defined here as the degree to which a system is likely to change due to exposure and sensitivity to a specific perturbation (Chapin et al. 2004). In this instance, vulnerability refers to the degree to which water resources, as represented by water budgets, and aquatic chemistry concentrations and loads, are likely to change due to exposure and sensitivity to permafrost thaw. Permafrost is ground (soil or rock and included ice and organic material) that remains at or below 0 °C for at least two consecutive years. It is defined by temperature rather than ice content and some permafrost, for example bedrock, may contain very little ice. However, permafrost commonly contains substantial amounts of ground ice. Permafrost thaw refers to the progressive loss of ground ice in permafrost, usually due to input of heat (IPCC 2019). Progressive thaw with active layer thickening, talik development, and rising of the permafrost base increases hydraulic conductivity and connections between ground and surface waters (e.g., Lamontagne-Hallé et al. 2018), and induces thermokarst where ice-rich ground thaws to produce depressions in the land surface through processes that include soil consolidation and settlement, slumping, gullying, and active layer detachment (Harris et al. 1988).

Environmental factors that influence the vulnerability of water resources to permafrost thaw can be categorized two ways: current landscape conditions and potential stressors (Fig. 1). Current conditions are a function of the landscape and its history and include permafrost and terrain characteristics. These control the sensitivity of water resources to permafrost thaw. When stressors are absent it is assumed that permafrost conditions remain static and there will be no permafrost thaw. Permafrost stressors include climate change and landscape disturbance. When these are present, permafrost warming and thawing, accompanied by thermokarst, occur over space and through time at varying rates.

Changes in permafrost caused by thaw and thermokarst are associated with several mechanisms that are associated, in turn, with changes in hydrologic and aquatic chemistry regimes. Studies documenting these changes provide insight into which influential factors (both conditions and stressors) were used in the CWRVI_PT (Table 1). Furthermore, interactions among hydrologic and aquatic chemistry regimes have cascading effects on aquatic ecology and ecosystem services, which are benefits humans directly gain from functions of the environment that provide, support, and regulate resources (Leemans and de Groot 2005). The interactions among these regimes also effect human water use for domestic and industrial purposes (Fig. 1). The focus of CWRVI_PT is on hydrological and aquatic chemistry regimes as defined by typical water budgets, aquatic chemistry concentrations, and loads.

The spatial domain selected for the development and application of the CWRVI_PT is the continuous and discontinuous permafrost region of Canada as defined using the circumpolar permafrost and ground ice data set of the National Snow and Ice Data Centre (http://nsidc.org/data/GGD318), a digitized version of the Circum-Arctic Map of Permafrost and Ground-Ice Conditions (Brown et al. 2002). Glaciers and large lakes (e.g., Great Bear Lake) were exempted because they are dominated by different processes than those represented using the methodological approach. Several large rivers (e.g., Mackenzie and Nelson) that pass through this domain have headwaters outside of the permafrost domain. Changes to streamflow and chemistry in such rivers that emerge from more southern regions may be due to a diversity of factors in addition to those accounted for by CWRVI_PT.

The temporal domain of the study is partially a function of data availability. Landscape and permafrost data were compiled over various periods of record with some dating back to the first decade of the 2000s. Climate change impacts are typically assessed using projected differences between a baseline or current period and a future period. In this study, these were 1986–2005 and 2030–2050, respectively. However, it is recognized that permafrost is already responding, and will continue to respond to ongoing and future warming and landscape disturbances. There will be spatial and temporal differences in permafrost and water resources responses to the expected change in climate between 1986–2005 (1990s) and 2030–2050 (2040s), with some impacts taking decades to emerge. Therefore, the exact permafrost and landscape conditions that will exist in the future are not known, and given this, future climate stressors were applied to current permafrost and landscape conditions. For these reasons, the temporal domain for the study was defined as from the 1990s to the 2040s, but it is recognized that both the earlier and later boundaries are uncertain.

The methodological approach to developing CWRVI_PT was to collate spatial data of indicators of landscape conditions and overlay these with spatial data of landscape disturbances as well as projected climate information. For each indicator that could be mapped across the domain, a ranking system was developed, which was then used to calculate vulnerability. Similar approaches were used by Allan et al. (2013) and Alessa et al. (2008). Following the framework described above, vulnerability was assessed as the sum of four sub-indices: permafrost conditions (I_p), terrain conditions (I_t), climate stressors (I_c), and disturbance stressors (I_d).

Each sub-index was a function of the ranks of several indicators. Each indicator was selected based on two main criteria: (i) that it has a significant impact on permafrost thaw and consequently water resources (Table 1), and (ii) that a pan-Canadian spatial data set exists (Table 2). The latter criterion eliminates some potential indicators, such as the influence of groundwater on permafrost thaw. Each indicator was selected because it has been documented as a factor controlling how water budgets or chemistry respond to permafrost thaw. How much control each indicator has on water resource vulnerability was ranked on a spectrum from relatively resilient (0) to vulnerable (1). Ranking divisions were also based on the classification, magnitude, or amount of variation of that indicator. Some vulnerability indices apply a weight to the ranks of each indicator and index based upon perceived influence (Allan et al. 2013). However, due to the difficulty of objectively determining which indices and indicators are more important, no weights were assigned to the ranks in this assessment. Instead, each sub-index was the average of the rank of each indicator so that no single indicator biased the overall index (see Alessa et al. 2008). When indicators and sub-indices had multiple components, the arithmetic mean of these components was calculated to ensure they were treated as a single combined value. This prevented biasing the index due to higher values associated with data availability.

The permafrost conditions sub-index was calculated following:

Ranks denoting permafrost extent (P_z) were associated with permafrost zones used in the Circum-Arctic Map of Permafrost and Ground-Ice Conditions (Heginbottom and Dubreuil 1993; Brown et al. 2002). This defines zones of continuous permafrost as those areas with 90%–100% spatial extent, discontinuous permafrost from 50% to 90%, sporadic from 10% to 50% and isolated between 0% and 10%. Although the concept of permafrost zones has been very useful at coarse scales, its major disadvantage is that it conceals finer scale heterogeneity. Furthermore, because of a lack of ground temperature measurements in many remote areas, zonal boundaries are also informed by air temperatures. These have increased since the late 20th century, introducing further uncertainty into the mapped zones. Gruber (2012) developed a permafrost zonation index from updated climate and terrain data at a finer scale of 1 km², which is applied here.

The potential of permafrost thaw to introduce new hydrological pathways and connectivity, or thermokarst features that change drainage networks (Table 1), is inherently related to its extent. In general, vulnerability due to permafrost thaw increases with decreasing permafrost extent (Zhang et al. 2008) because the permafrost becomes warmer and thinner and it is likely to cause more landscape changes than merely an increase in active-layer thickness (as occurs in continuous permafrost). This spatial dimension has implications for water flow pathways, streamflow seasonality and chemistry (Petrone et al. 2006; O’Donnell et al. 2012; Connon et al. 2014) and has been attributed to differences in rates of lake drainage (Smith et al. 2005; Karlsson et al. 2012). This implies that where the permafrost zonation index is between 50% and 90% (i.e., discontinuous) water resources are most vulnerable to permafrost thaw. The next most vulnerable areas were identified as those where the index ranges from 10% to 50% (i.e., sporadic). Although thaw may not result in complete disappearance of continuous permafrost, even subtle changes in active-layer depth can have profound impacts on water chemistry (Lewis et al. 2012; Lafrenière and Lamoureux 2013; Farquarson et al. 2019). This can have widespread effects on water resources because the frozen ground is so extensive (>90%) (Tables 3 and 4). Where there is almost no likelihood for permafrost (<10%, i.e., isolated), we assumed that its thaw will not be significant to water resources response at the landscape scale.

Ground ice content (C_i) ranks were based on estimates of relative abundance of excess ground ice from O’Neill et al. (2019) (Table 2). Ground ice content controls the amount of water that is hypothesized to be a source of increased streamflow by providing moisture inputs that were previously bound as ice (Connon et al. 2014). These inputs are expected to manifest as higher subsurface flows. The presence of wedge ice increases the likelihood of lateral breaching of lakeshores, catastrophic lake drainage (Pohl et al. 2009; Jones et al. 2011), and subsequent changes in surface storage capacity and drainage network structure. These findings imply those areas with higher ground ice content have water resources that are more vulnerable to permafrost thaw (Table 3). This is a similar approach to that used by Olefeldt et al. (2016) in a thermokarst landscape classification study. There are three common types of excess ground ice: relict, segregated, and wedge, and each can occur in high, medium, low, and negligible fractions (O’Neill et al. 2019). Ice abundance is not equal for each fractional category. For example, areas of high density wedge ice abundance may occupy 3.5%–16% of the upper 5–6 m of permafrost (Pollard and French 1980; Couture and Pollard 1998), whereas segregated ice may commonly exceed this volume (Pollard and French 1980). Relict ice, where present, may represent areas of highest relative ground ice abundance (Couture and Pollard 1998) as evidenced by extensive retrogressive thaw slumping (Kokelj et al. 2017) owing to its high ice content and spatial distribution. Vulnerability was assessed for the relative contribution of each ground ice type according to its abundance (high, medium, low, and negligible) (Table 4) whereby the maximum value of the summed contributions was equal to 1. The actual ranks were summed for each type of ground ice so that C_i is equal to:

where W_i, S_i, and R_i are the wedge ice, segregated ice, and relict ice ranks, respectively.

The terrain sub-index was calculated as:

Data to assess overburden (O) and surficial geology (G_s) ranks came from Fulton (1995) (Table 2), which details the mapping of surficial geology across Canada. Two classes of overburden thickness were recognized: (i) thin (veneer), and (ii) thick (all other types). Thin overburden indicates a lower potential for significant landscape and water resource regime changes with permafrost thaw (Lamoureux and Lafrenière 2009; Dugan et al. 2012; Kokelj and Jorgenson 2013), and thus was considered resilient (ranked as 0; Table 3). Thick overburden was ranked as 1 (Table 4). The influence of surficial geology on water resources vulnerability is a function of hydraulic conductivity and particle size (Smith and Burgess 1999; Olefeldt et al. 2016) as subsurface flow can accelerate permafrost thaw (McKenzie and Voss 2013; Lamontagne-Hallé et al. 2018). For example, offshore glaciolacustrine deposits are finer than glaciofluvial deposits, with the former provided a ranking of 0.2 and the latter 1. Sands and gravels from weathered bedrock were ranked as more vulnerable (0.7) than low conductivity till deposits (0.5) (Table 4). Similarly, higher conductivity associated with bedrock geology (G_b) would result in a larger response in subsurface pathways, rates of subsurface flow, and potential changes in chemistry as pore space and fractures become available to transmit water with permafrost thaw (Carey and Woo 2001; Walvoord and Striegl 2007). Using Freeze and Cherry (1979) and bedrock geology data from Wheeler et al. (1997) higher hydraulic conductivity bedrock was deemed relatively more vulnerable (Table 4).

Water resources in areas with thicker organic layers are expected to be more vulnerable to change because the thaw of any frozen ground in these areas would result in more profound changes to water budgets and chemistry than in areas with thinner organic layers (Carey and Woo 2001; Buffam et al. 2007; Walvoord and Striegl 2007) (Tables 3 and 4). As no pan-Canadian database of organic layer thickness exists, soil organic carbon content from the Northern Circumpolar Soil Carbon Database (Hugelius et al. 2013) can be used as a proxy because it should correlate with thicker organic layers (Olefeldt et al. 2016). The required soil organic carbon content ranks (C_s) were calculated from bulk density, organic carbon percentage, soil depth, and fraction of coarse fragments and (or) segregated ice content characteristics of representative soil types in northern Canada and extrapolated by appending data to an existing soil type map of Canada (Tarnocai and Lacelle 1996).

Assessing vegetation-permafrost-hydrology-biogeochemistry interactions must include a representation of land cover (Jorgenson et al. 2013). Evaporation rates can differ among land cover types associated with permafrost extent and ground ice type (Wright et al. 2009; Raz-Yassef et al. 2017). If vegetation type distribution associated with permafrost thaw is extensive enough, regional evapotranspiration rates could be altered. The role of land cover on water resource vulnerability to permafrost thaw was assessed using a product from the Canada Centre for Remote Sensing (CCRS) (Table 2). The original 31 categories in the database were re-classified into four: (i) forests, (ii) tree line transition, (iii) shrubs and low vegetation cover, and (iv) exposed rock. Wetland presence was calculated using the Global Surface Water data set (Table 2) because although the CCRS data included several wetland categories, these include both wetlands and shrubs and so were not considered to be appropriate for northern Canada. The land cover rank, L, was calculated using ranks of each land cover. Forests were ranked as less vulnerable to permafrost thaw-induced water resource change as they protect permafrost by shading the ground surface (Jorgenson et al. 2010). The tree line transition was identified as a separate land cover type because of its distinct mixture of tundra, trees and, importantly, expanding shrubs that can influence ground thermal conditions (Liston et al. 2002; Lawrence and Swenson 2011). Exposed bedrock was ranked the least vulnerable of the land cover types. Wetlands can be particularly vulnerable (Wright et al. 2009; Quinton et al. 2011; Connon et al. 2014) so wherever the Global Surface Water data set indicated wetland presence, another 0.3 was added to L (Tables 3 and 4).

The only portion of Canada that contains yedoma (organic rich permafrost with 50%–90% ice) is in central Yukon (Wolfe et al. 2009) (Table 2). Areas with yedoma (Y) were ranked highly vulnerable (i.e., 1) because of higher likelihood of excess ground ice (Tables 3 and 4) and extensive thermokarst, and thus a high potential for significant landscape change (Strauss et al. 2017). The large contiguity of yedoma regions in the Yukon make them significant due to the potential for widespread re-mobilization of organic carbon (Strauss et al. 2017), which can impact regional-scale water quality when permafrost degrades.

Relief (R) categories and data were obtained from Gruber (2012) (Table 2). The most vulnerable areas (i.e., 1) were flat (zonation index of 1.5–0) (Tables 3 and 4) because of possible subsidence, inundation and changing drainage network structure. This may also alter subsurface pathways with implications for water chemistry (Sannel and Kuhry 2011). Hilly and undulating areas (zonation index from 3.5 to 1.5) were intermediately ranked as they are more susceptible to hillslope thermokarst processes and associated changes these may have on water chemistry. Rugged areas were ranked less vulnerable (Table 4) as they are less likely to have significant amounts of ground ice (Gruber and Haeberli 2009), as are wetlands or lakes (Olefeldt et al. 2016).

This sub-index (I_c) focuses on stress due to projected climate change on permafrost and associated impacts on water resources. The climate data consists of monthly averages of precipitation and temperature from an ensemble of 29 Coupled Model Intercomparison Project (CMIP)5 climate models using the representative concentration pathway (RCP) 8.5 emission scenario, which represents a “business as usual” scenario based on current global greenhouse gas emissions (Environment and Climate Change Canada 2018). The data sets are at a 1° spatial resolution across Canada and matched to the land areas, but are bounded between 50° and 80°N and 50° and 140°W, which reduces the size of the domain from the maximum extent of permafrost on Canada’s land mass. The 50th percentile of the 29 CMIP5 models was employed in our assessment to represent the average projection from all climate models. Individual indicators were calculated from monthly, seasonal, or annual precipitation and temperature data. All climate ranks were a function of the change between the near future (2031–2050) and current conditions (1986–2005). For the purposes of sub-dividing the indicators and providing rank values, the range of climate indicators was based on the greatest amount of change between current and near-future values. For each indicator, 11 quantiles were calculated, using the first and last to isolate the upper and lower extremes (Table 2). The range between the two quantiles was equally divided into 9 ranges to fill the rankings from 0.1 to 0.9 (Table 4). I_c was calculated for the near future period as:

where:

and

Air temperature (T) and precipitation (P) were the two selected climate parameters, with specific indicators within each. Four key aspects of temperature change were considered, which affect the amount and timing of heat exchange between the atmosphere and the land surface (Table 1). An increase in mean annual air temperature (MAAT) results in initiation of permafrost thaw, increasing summer thaw depth, deepening of the active layer, and reduction of the extent of frozen ground (Rinke et al. 2012). Warmer winter temperatures (average of December to February; T_w) diminish ground cooling necessary to maintain permafrost and can prevent complete freezing of the active layer, which is especially important over warm permafrost within which connectivity can be maintained (Smith et al. 2012; Connon et al. 2014). An increase in summer temperature (average of June to August; T_s) can be associated with higher evapotranspiration rates. It also increases ground temperatures and active-layer thickness and with all else being equal will lead to higher ground temperatures the following winter (Rinke et al. 2012). Warmer soil is associated with enhanced chemical cycling (Koch et al. 2013) and a higher likelihood of subsurface hydrological connectivity (Walvoord and Striegl 2007). Extremely warm summer temperatures can precondition areas for, or trigger, thermokarst processes, such as active-layer detachments that impact water chemistry (Dugan et al. 2012; Lewkowicz and Way 2019). Longer open-water season (annual ice-free duration of waterbody surfaces; D) permit higher heat storage in lakes, which may delay winter freeze-up dates (Barnhart et al. 2016), permitting higher rates of heat conduction between waterbodies and the adjacent ground. This heat conduction can increase near-surface ground temperatures in areas below water bodies by up to 10 °C more than MAAT (Jorgenson et al. 2010). These examples demonstrate how multiple temperature indicators are necessary to capture the different ways in which climate warming can influence permafrost thaw and subsequent water resource response. Ranks for MAAT, T_w, T_s, and D were all calculated similarly. The greatest increase was considered the most vulnerable (1) and lowest change the least vulnerable (0), with sub-rankings within this range calculated following the quantile method described above. Ranks for D were established in the same way with values approximated by determining the duration of each calendar year when air temperatures were above 0 °C.

Precipitation is an important stressor on permafrost as it generates potential for soil saturation, low basal shear strength, and erosion (Lewkowicz 2007; Lamoureux and Lafrenière 2009; Kokelj et al. 2013), or in the form of snow, it acts as an important ground insulator (e.g., Stieglitz et al. 2003). The two indicators included were late summer rain (August precipitation; P_A) and annual snowfall (P_W). Rain is more likely to influence runoff chemical concentrations and loads during late summer when thaw depths are greatest. This increases the likelihood of rapid permafrost thaw (Rinke et al. 2012), rapid drainage of thermokarst lakes (Pohl et al. 2009) or hillslope thermokarst formations (Lewis et al. 2005; Lewkowicz and Harris 2005; Lacelle et al. 2010), each of which can have significant sudden impacts on water budgets, sediment concentrations and loads to nearby waterbodies (Lewis et al. 2012; Lafrenière and Lamoureux 2013). These are responses to individual runoff events, when rain is heavy and evaporation is limited, but changes in drainage networks or movements of large volumes of sediment can have long lasting repercussions for water chemistry and streamflow production (Kokelj et al. 2013). August precipitation was selected because it is the latest month in which precipitation could be reasonably assumed to be entirely composed of rainfall across the entire spatial domain. P_A was ranked such that an increase indicates permafrost is more vulnerable, whereas a decrease may indicate more resilience to change. The greatest increase in rainfall was ranked most vulnerable (1) and any decrease or no change in rainfall ranked as least vulnerable (0) with sub-rankings determined by the same quantile approach as with temperature.

The length of the period during which air temperature is below 0 °C was used to estimate when precipitation would fall as snow, a necessary step because the climate data set does not differentiate precipitation phase. Annual approximate snowfall (P_W) calculated in this way was summed and averaged for each time period. The quantile approach was applied, but because there are both projected increases and decreases in P_w, quantiles were aligned so that no change was given a value of zero. Snowfall is an indicator of ground insulation during winter (Rinke et al. 2012) with higher amounts reducing ground cooling (Stieglitz et al. 2003) leading to increasing thaw depths in spring and summer (Walvoord and Kurylyk 2016), which has implications for opening and sustaining subsurface pathways for runoff. Less snow is expected to increase the resiliency of the permafrost. Proportional increases and decreases in snowfall were used with the greatest increase at 0.7, the greatest decrease at −0.3, and no change at 0.

The disturbance sub-index was calculated following:

where HP is human presence (HP = 1 if presence = “yes”, otherwise = 0) and FP is fire potential as:

The HP disturbances considered were roads, mines, oil and gas wells, pipelines, communities, and distant early warning (DEW) line sites (Table 2). The response of water resources, specifically runoff chemistry, to these disruptions depends on the degree of permafrost thaw, particularly in discontinuous permafrost (Williams et al. 2013). Linear disturbances, such as roads (Huntington et al. 2007; Wolfe et al. 2015), pipelines (Smith et al. 2008), and seismic lines (Williams et al. 2013), affect the ecosystem that may be providing protection for permafrost. Disruptions include vegetation clearing, overburden removal, drainage alteration, artificially warming or cooling the ground, and surface water management. These activities melt ground ice, change snow accumulation, and alter drainage networks through subsidence. For example, soil moisture commonly increases along roads due to blockage of runoff by embankments and the energy advected with this water results in earlier active layer thaw (Williams et al. 2013). Mines and oil and gas wells are areas of increased permafrost vulnerability, as vegetation is cleared, and overburden is commonly removed allowing for more direct radiative fluxes that encourage permafrost thaw (Cao et al. 2016). Active on-site water management for operations changes water surface pathways, chemistry, and runoff rates. Pipeline corridors are another linear disturbance as permafrost thaw and settlement occurs due to alteration of the ground thermal regime (Smith et al. 2008), with implications for water resources, specifically water budgets. The impact of constructing a right-of-way for a pipeline or seismic line on ground temperatures can outweigh climate warming impacts for several years after construction (Smith and Riseborough 2010). Communities impact permafrost through vegetation changes, surface disturbances, and development (Huntington et al. 2007), whereas DEW line sites may also cause contamination (Thomas et al. 1992; Poland et al. 2001). Disturbed areas were ranked as vulnerable (i.e., 1) with all other areas ranked as 0. Disturbed locations were imported as point or line data and subsequently buffered by 100 m. No attempt was made to predict disturbances associated with future development.

The influence of climate warming will dominate that of forest fire to reduce permafrost extent over time (Zhang et al. 2015), but forest fire can reduce thin, discontinuous permafrost extent by 16% or more (Yoshikawa et al. 2002; Gibson et al. 2018). Fire accelerates permafrost thaw because the loss of insulation with organic soils allows for a thicker active layer that is less likely to recover as the climate warms (Yoshikawa et al. 2002; Holloway et al. 2020). The deepening of the active layer, or loss of permafrost altogether with fire, opens previously unavailable pathways to groundwater that can change streamflow seasonality and chemistry. Forest fire disturbance was represented in CWRVI_PT with a form of fire weather index (Van Wagner 1987) based on summer temperature (T_s) and summer precipitation (P_s). Its application was limited spatially to forested areas. Using the same climate data described above for use in I_c, the change in summer precipitation between time periods was compared, and areas with the smallest increase were ranked as most vulnerable (P_s = 1) and areas with greatest increase were most resilient (P_s = 0). The same method to derive values of summer temperature ranks (T_s) using quantiles used from the I_c for summer temperature was used for FP. The arithmetic mean of the two indicators constituted FP (eq. (9)).

Figure 2a illustrates the spatial distribution of the values of permafrost conditions (I_p), which increase northward from the southernmost edge of permafrost. The highest values occur across a wide swath extending from Ungava Bay in Québec through the Hudson Bay Lowlands, northwest past the southern shores of Great Bear Lake towards the Peel Plateau in Northwest Territories. A portion of the Mackenzie Mountains has relatively lower values because of low ground ice. Variation within these predominant spatial trends is due to the presence of water bodies and ground ice. Regions with continuous permafrost and low ground ice (e.g., large parts of Baffin Island) contrast regions with extensive areas that contain all three types of massive ground ice (e.g., Banks and Victoria Islands) (O’Neill et al. 2019), which enhances values of CWRVI_PT. This highlights the result that the presence of continuous permafrost does not necessarily result in lower values of CWRVI_PT.

Calculation of terrain conditions (I_t) incorporates the greatest number of indicators of all the sub-indices, which results in the most spatial variation although broad scale patterns exist (Fig. 2b). Values of I_t are lowest in the mountainous regions of the three territories and Labrador because of high relief and widespread exposed bedrock. The moderate relief, thin overburden, and exposed bedrock landscapes of the Canadian Shield in mainland Nunavut and the Ungava Peninsula in northern Québec exhibit values from 0.4 to 0.5. Portions of Baffin Island display values as high as 0.6 due to thicker overburden and low relief. This combination of characteristics that result in values that can exceed 0.7 also occurs across much of the Hudson Bay Lowlands and the Mackenzie River Valley in the Northwest Territories.

The sum of I_t and I_p, which could be considered as the sensitivity of water resources to permafrost thaw, is highest in two regions that have values approaching 1.5. The first extends from the northern Yukon south towards the southern Northwest Territories. The second is in northern Manitoba. The highest values generally occur in flatter regions of warm ice-rich discontinuous permafrost with widespread organic soils. Values around 1.0, in northern regions of Québec and Labrador are comparable with those in northern Alberta and Manitoba. Values in the Ungava Peninsula and much of mainland Nunavut are lower because the thinner coarser overburden contains less ground ice and is less susceptible to thaw-induced thermokarst (Olefeldt et al. 2016). Combined values of I_p and I_t as low as 0.2 occur in the mountainous regions of northern Baffin Island. In contrast, Banks and Victoria Islands, at a similar latitude, have higher values due to high ground ice contents.

The influence of climate stress on water resource vulnerability to permafrost thaw is most pronounced on Victoria Island, northern Baffin Island, and the region south of the Queen Maud Gulf. In these areas, I_c exceeds 0.6 (Fig. 2c) primarily because of higher winter temperatures and greater precipitation. I_c values tend to decrease towards Ellesmere Island, and to the southwest towards the Northwest Territories and Yukon. Instances that deviate from this trend occur in northern Québec, the Hudson Bay coastline of Manitoba, northwest Alberta, areas west of Great Bear Lake, and the northern Yukon. These areas exhibit high I_c values that range from 0.5 to 0.7.

Stressors associated with HP can only be detected at the scale of the entire study domain in a few locations. One of these is in northern Alberta where the network of pipelines and roads is extensive (Fig. 2d) and values of I_d are as high 0.7. It is possible that the ecosystem processes that protected permafrost in northern Alberta have already been removed from these locations of landscape disturbance, so the impacts on water resources may have already happened. The FP is zero in arctic and alpine landscapes, so stressors identified in these areas are limited to locations with HP disturbances. More frequent forest fire disturbance is predicted south of 60°N particularly in northern Saskatchewan and Manitoba, and in Ontario west of James Bay, where FP is at a maximum of 0.4.

Figure 3 illustrates the results of applying eq. (1) to assess the vulnerability of water resources via permafrost thaw into the middle of the 21st century (from where data can be downloaded is described in Supplementary File S1¹). The patterns exhibited in Fig. 2 are projected to continue into the future with enhanced vulnerability across a wide swath of the Northwest Territories, and northern Manitoba, and Ontario. Values of CWRVI_PT indicate regions of higher vulnerability south of Hudson Bay, portions of central Nunavut, and southwest of Great Bear Lake (Fig. 3). Relatively high CWRVI_PT values in central Baffin Island, western Victoria Island, and portions of central Nunavut can be attributed to high I_c (Fig. 2c). The Mackenzie Mountains, and portions of the Boothia Peninsula are relatively less vulnerable. Mountainous regions of Baffin Island also have low CWRVI_PT values, due to very low values of I_t, I_p, and I_d. In contrast, the I_c values in these regions are among the highest (0.75) and dominate final CWRVI_PT values (see Fig. 2c). The high I_c values are due to projected increases in winter air temperatures and summer precipitation. Vulnerability values averaged 0.8 ± 0.3 using only a combination of I_t and I_p, but the average increased to 1.3 ± 0.3 when climate and disturbance stressors were added. The increase across the entire pan-Canadian domain was almost exclusively attributable to climate stress.

The values assigned in Table 4 are admittedly only informed estimates. For this reason, an evaluation of CWRVI_PT is important. However, this is challenging because there is a paucity of large-scale data concerning water resource change across Canada that could be used to evaluate CWRVI_PT at scales beyond individual headwater catchments or landscapes. Studies of this type have taken place in Siberia and Alaska using satellite observations or widespread water chemistry monitoring networks (e.g., Smith et al. 2005; Lyon and Destouni 2010). The few comparable studies in Canada that discuss implications on water resources include Li Yung Lung et al. (2018) who, in a synoptic analysis of chemical loads from the Canadian land mass to the Arctic Ocean, demonstrated that sparse regional aquatic chemistry data make evaluation difficult. Lantz and Kokelj (2008) detected an increase in thaw slump density across the uplands adjacent to the Mackenzie Delta in the latter half of the 20th century. Similarly, Marsh et al. (2009) documented changes in the rate of catastrophic lake drainage associated with thawing and thermokarst across an area east of the Mackenzie Delta. Follow-up research documented widespread water chemistry impacts (Deison et al. 2012) but these results are geographically limited to one landscape among many within the study domain. Lantz and Turner (2015) determined that thermokarst and catastrophic drainage were responsible for half the lakes that experienced a change in area in a portion of Old Crow Flats, Yukon between 1951 and 2007. They suggested that further research was needed on the mechanisms and conditions that cause catastrophic drainage to extrapolate findings wider afield. This also demonstrates that evaluation data are limited.

One approach to systematically evaluate CWRVI_PT is to select the most sensitive areas based on permafrost and terrain conditions and examine if these have experienced changes in water resources. The most sensitive areas are those in the top quartile of the combined value of I_t and I_p (i.e., >0.964). Where the sum of I_t and I_p exceeded 0.964, the higher of the two sub-indices was mapped (Fig. 4). Permafrost conditions generally matter more in the discontinuous permafrost zone with areas of medium ground ice abundance than terrain in the continuous permafrost zone where ground ice abundance is high. Expected water resource responses in these regions include changes in drainage network structure and surface storage due to ground ice melt and subsidence as well as changes in chemical fluxes associated with mass movement (Table 1). These responses are especially associated with thermokarst development. Areas where I_p dominates CWRVI_PT align well with areas of high retrogressive thaw slump density and changes in water chemistry in the northern Yukon, northern mainland Northwest Territories, and portions of Banks and Victoria Islands (Kokelj et al. 2017). These locations are associated with the margins of the late Wisconsinan ice sheet and have high concentrations of relict ice (Kokelj et al. 2013). Thaw slumps increase total suspended solids by orders of magnitude and double specific conductivity in water bodies immediately downslope, and can increase summer SO₄/Cl ratios hundreds of kilometers downstream. Such thaw slumps also have the potential to cause substantial change to the biological composition and food web structure of lakes and rivers (Moquin et al. 2014; Chin et al. 2016).

Terrain plays a dominant role in water resource vulnerability where relief is limited, wetlands are common, and there is high soil organic carbon content (Fig. 4). These factors control the active layer thickness and the areal extent of frozen ground, altering runoff pathways by opening previously inaccessible subsurface routes. Runoff chemistry changes as minerals and solutes available to runoff are impacted (Table 1). Scotty Creek watershed, southeast of Fort Simpson, Northwest Territories, is within a region in which CWRVI_PT values are near their highest and I_t values dominate (Fig. 4). This is an area of thick organic soils, flat terrain, numerous wetlands, and discontinuous permafrost with specific land cover types associated with permafrost presence. Bogs and fens tend to be permafrost-free, whereas peat plateaus contain permafrost, and the latter decreased 40% between 1947 and 2008 throughout the Scotty Creek basin (Quinton et al. 2011). Loss of this land cover type and the permafrost associated with it can open up previously inaccessible portions of the subsurface to groundwater flow (Connon et al. 2015). Associated with this extensive permafrost thaw, and changing structural hydrological connectivity, regional runoff ratios and annual streamflow have been increasing steadily since 1970 at a rate between 4% and 10% annually due to landscape changes and not precipitation (Connon et al. 2014).

A second example of links between I_t and water resource vulnerability concerns permafrost thaw and runoff in winter. St. Jacques and Sauchyn (2009) hypothesized that permafrost thaw was responsible for widespread increases in winter baseflow across the Northwest Territories, citing work by Smith et al. (2007) and Walvoord and Striegl (2007) that implied that this shift must be due to higher flows along better connected groundwater pathways. Limited data show a strong positive correlation between I_t and changes in winter baseflow between the mid-1970s and mid-2000s (r = 0.91) (Fig. 5).

I_c and I_d were not as extensively predominant in the CWRVI_PT as I_t and I_p. However, I_c was the predominant sub-index in portions of the High Arctic and northern Québec (see Fig. 2c). These areas could be interpreted as more resilient to climate change because greater change in climate is required to induce vulnerability due to colder and more extensive permafrost. In areas where permafrost and terrain conditions are the predominant controls, less climate change may be required to induce an impact or exhibit tipping points, and they may already be exhibiting widespread changes in water resources (e.g., St. Jacques and Sauchyn 2009; Kokelj et al. 2017). I_d was the highest sub-index only along the southern reaches of the permafrost zone in locations of HP (see Fig. 2d). This indicates changes to water resources associated with permafrost thaw in these areas may be dependent on land use decisions (e.g., seismic line intensity).

Some potentially important factors could not be included or were represented by proxies (e.g., soil organic carbon content) in CWRVI_PT due to a lack of data. For example, groundwater discharge areas are more prone to permafrost thaw and are especially important for maintaining water levels within wetlands and lakes (Yoshikawa and Hinzman 2003). These are not represented in CWRVI_PT, which therefore underestimates vulnerability in these areas. Seismic cut lines created during hydrocarbon exploration are also known to lead to permafrost thaw (Williams et al. 2013), but were excluded to prevent spatial bias linked to incomplete data coverage.

Figure 3 represents a first attempt at characterizing the vulnerability of water resources to permafrost thaw across northern Canada. The caveats described above, especially the basic information gaps, must be considered when interpreting and applying the index. Further evaluation is needed at a variety of locations across a spectrum of CWRVI_PT values to confirm if higher values are associated with greater impact or larger changes in water resources. Larger changes should be reflected in metrics, such as stronger trends or step changes in streamflow or solute concentrations. The CWRVI_PT could be developed further in an iterative manner. The first step would be to identify expected hot spots of change and (or) uncertainty from existing results and prioritize those areas for observation and research. The results of this field research could inform and improve ranking schemes used to derive CWRVI_PT. A second step would be to apply numerical hydrological and biogeochemistry models to see if they reproduce changes that align with spatial patterns of CWRVI_PT. Risk results from the interaction of vulnerability, exposure, and hazard (Oppenheimer et al. 2014) and is often expressed as the probability of occurrence of hazardous events, or trends multiplied by impacts if these events occur. Therefore, CWRVI_PT could represent a measure of vulnerability used in an analysis of specific risks to communities, aquatic species or ecosystems, environmental services, infrastructure and natural resource developments in northern Canada.