Active layer variability and change in the Mackenzie Valley, Northwest Territories between 1991-2014: An ecoregional assessment

ABSTRACT

Active layer thicknesses (ALTs) from sites along a transect through the Mackenzie Valley, Northwest Territories, Canada, were analyzed to explore variation in thickness within and between ecoregions. At an ecoregional scale the relation between ALT, latitude, freezing and thawing degree-days, and snowfall were examined to determine the presence of trends. Site-specific variables including dominant vegetation and substrate were explored to explain spatial variability in ALT within ecoregions. Generally, average ALT increases moving southward through the comprising ecoregions (68 cm to 126 cm), following the increase in air temperature. Spatial variability in ALT within ecoregions was greater than that between ecoregions (up to 145 cm), which may be attributed to site-specific conditions (vegetation and snow cover). Most notable, sites with shrubs had thicker than average active layers likely because of increased snow retention leading to warmer overall ground conditions. Despite a warming trend in air temperatures, only one northern ecoregion showed a corresponding thickening trend in ALT. Sites located in southern ecoregions with mature forests showed limited response to changes in air temperature. For these locations, disturbance, specifically changes in thermally protective vegetation cover, rather than changing air temperature could potentially have a larger impact on ALT into the future.

KEYWORDS:

• permafrost
• active layer
• ecoregion
• Mackenzie Valley

Introduction

The active layer refers to the layer of ground in areas underlain by permafrost, which is subject to annual freezing and thawing and is the location where most subsurface ecosystem processes occur in the permafrost landscape. The nature of this layer is dynamic, responding to temporal climatic shifts showing only subtle connection to atmospheric interannual variability (Shur, Hinkel, and Nelson 2005). Although the most common metric used to describe the nature of the active layer is thickness at the time of maximum thaw, the predictability of this value is not universally tied to latitude or even mean annual air temperature (Bonnaventure and Lamoureux 2013). Active layer thicknesses (ALTs) are dependent on several factors, including latitude, elevation, solar insolation patterns, albedo, proximity to glaciers, and/or warm maritime ocean currents, as well as site-specific factors including soil texture and moisture content (Bockheim and Hall 2002). Some factors including vegetation and snow distribution can vary substantially within climatic zones due to topography, hydrology, and dominant wind patterns (Evans et al. 1989; Young, Woo, and Edlund 1997; Garibaldi, Bonnaventure, and Lamoureux 2021). Vegetation has been shown to influence ALT through snow entrapment and interception, shading, evapotranspiration, and insulation depending on the vegetation type (Duchesne et al. 2015; Fisher et al. 2016; Loranty et al. 2018). Previous studies of ALT have primarily considered the impact of summer (thawing degree-days) conditions on ALT; however, antecedent winter conditions including snow cover and freezing degree-days may also substantially impact ALT (Mazhitova et al. 2004). Seasonal snow cover insulates the surface from cold winter conditions, leading to warmer winter ground conditions and potentially thicker active layers (Zhang 2005). Snow cover also alters the timing of freeze-back in the autumn, the progression of the ground freezing front, and the nature of the spring thaw (Zhang 2005; Morse, Burn, and Kokelj 2012). Thicker snow covers can also have considerable impact on the ground surface temperature by trapping latent heat expelled from the subsurface that cannot be freely dispersed to the atmosphere (Zhang 2005; Bonnaventure, Lewkowicz, and Hillaire-Marcel 2008).

The complexities of the active layer are further complicated by conflicting definitions of the active layer itself (Burn and Smith 1988). Thermally, the base of this layer ideally exists at the idealized freezing/melting point of water in nature (0°C); however, this may not match the base of the physical active layer. Much of this discrepancy is attributed to the intricacy and diversity of active layer materials, which can contain variable levels of moisture and substrate composition including mineral sediment, ions in solution, and organic layers (Hinkel et al. 2001; Shiklomanov et al. 2010). Additionally, the active layer is morphologically linked to the environment and can be influenced by the adjacency to water bodies or mass movement debris (Hinkel and Nelson 2003).

Because the active layer can be spatially heterogeneous on multiple scales, assessing variability is challenging but necessary to better understand the impact of change in cold regions (Bonnaventure, Lamoureux, and Favaro 2017; Garibaldi, Bonnaventure, and Lamoureux 2021). Understanding this variability helps to recognize the range and complexity that goes with trying to ascertain how permafrost environments will be impacted by warming (Bonnaventure and Lamoureux 2013). The goal of this study is to examine the temporal and spatial variability in ALT using an ecoregional approach. Recent studies have focused on the impact of vegetation and ecosystem characteristics on ALT largely focusing on tundra and boreal forests (Duchesne et al. 2015; Fisher et al. 2016; Loranty et al. 2018). Expanding these regional ALT assessments specifically in the Mackenzie Valley on an ecoregional scale will aid in small-scale (circumpolar, national scale) modeling studies which utilize regional data (e.g., Obu et al. 2019) and may allow for the creation of ALT models in this region (e.g., Yi et al. 2018). This study updates previous work in the Mackenzie Valley relating ALT to ecoclimatic zone and further analyzes an ALT data set used in previous work (Tarnocai, Mark Nixon, and Kutny 2004; Duchesne et al. 2015) building to explore the relation between active layer thickness, broad climate, and vegetation characteristics. The different ecoregions in the Mackenzie Valley have shown variable attributes with respect to the nature and timing of ground freeze-back in the autumn (S. L. Smith et al. 2016). Therefore, it is possible that active layer thickness also varies within and between these ecoregions. In order to conduct this analysis, multiple years of data must be examined and compared. This study utilizes data collected since 1991 from an active layer monitoring network in the Mackenzie Valley and Delta region to document ALT through time. Sites are representative of the range of vegetation and surficial material conditions throughout the region. The Mackenzie Valley is well suited for this type of study because both permafrost and vegetation are heterogeneous across the study transect, which passes through discontinuous (underlying <90 percent) to continuous (underlying 90–100 percent) permafrost zones and from boreal forested to tundra vegetation ecoregions. This study compares regional and temporal ALT patterns between ecoregions spanning the valley, expanding on work previously conducted by Duchesne et al. (2015) by further examining these data and conditions according to ecoregions.

Study area

The Geological Survey of Canada active layer monitoring transect spans the length of the Mackenzie Valley and Delta, with sites from the Delta and Tuktoyaktuk Peninsula to the southern Northwest Territories (Figure 1). Sites are in Inuvialuit, Gwich’in, Sahtu, and DehCho First Nations settlement regions and traditional territories. The transect (Figure 1 and Table 1) comprises multiple monitoring sites and extends from the zone of discontinuous to continuous permafrost and from the boreal forest near Fort Simpson to the Beaufort Sea coast (S. L. Smith, Riseborough et al. 2009). Substrate along the transect is variable but consists of mainly moraine, lacustrine, fluvial, and deltaic deposits that can be ice-rich, although ice content varies from site to site (Heginbottom, Dubreuil, and Harker 1995; Aylsworth, Burgess et al. 2000; Aylsworth, Kettles et al. 2000). Extensive peatlands and fens consisting of sedges, mosses, and stunted black spruce are particularly prevalent in the Southern portion of the valley (Aylsworth, Burgess et al. 2000; Aylsworth, Kettles et al. 2000; Ecosystem Classification Group 2007). The sites in the transect are representative of the various ecoregions in Northwest Territories, which are generally defined based on unique regional climatic gradients and major changes in physiography (Ecosystem Classification Group 2010). Ecoregions are classified furthermore in (Level II) ecoregional classification. This Level II classification adds an additional level of complexity based on macroclimate, vegetation, soil, geology, and physiographic features, which is the subdivision basis of permafrost and vegetation characteristics and additionally combined with climatic parameters to create the Level III ecosystem classifications (Ecosystem Classification Group 2012). Level III ecoregions that contain active layer monitoring sites in the Mackenzie Valley along this transect include Tundra Plains Low Arctic north (LAn), Taiga Plains High Subarctic (HS), Low Subarctic (LS), and Mid Boreal (MB), Boreal Cordillera High Boreal (HB), and Taiga Cordillera Low Subarctic (LS) encompassing seven of the nineteen possible Level III ecoregions found in the Northwest Territories (Figure 2 and Table 1). Throughout the valley, air temperature, annual precipitation, and snowfall increase moving north to south (Table 2).

Table 1. Summary table of characteristics of each ecoregion with active layer monitoring sites (Ecoregions Working Group 1989; Ecosystem Classification Group 2007, 2010, 2012; Environment Canada 2017).

Ecoregion	Latitude (°N)	Longitude (°W)	AMAT (°C)	Precipitation (mm)	Characteristic vegetation	Permafrost condition	Community
Tundra Plains LAn	66–70	118–130	−11	130–190 50% falling as snow	Low shrubs, sedge, and mosses	Continuous	Tuktoyaktuk
Taiga Plains HS	65–69	118–135	−5 to −11	250–340 50% falling as snow	White spruce (Picea glauca) (Krummholz)	Continuous and extensive discontinuous	Inuvik
Taiga Plains LS	62–66	117.9–132.8	−3.5 to −9	230–350 50% falling as snow	Open canopy white (Picea glauca) and black spruce forests (Picea mariana)	Extensive discontinuous	Norman Wells Fort Good Hope
Taiga Plains MB	60–62	112–123.5	−2 to −5.5	310–410 45% falling as snow	Mixed wood forests of aspen (Populous tremuloides) and white spruce	Extensive to sporadic discontinuous	Fort Simpson
Boreal Cordillera HB	61–62	123–126	−4 to −5	340–400 40% falling as snow	Lodge pole pine (Pinus contorta)	Extensive discontinuous	Wrigley
Taiga Cordillera LS	62–64	123.7–130	−7	280–350 40% falling as snow	Open spruce woodlands and Jack Pine (Pinus banksiana)	Continuous in the northern half of the ecoregion to extensive discontinuous in the southern portion	None

Table 2. Climate Normals (1981–2010) for five of the Environment Canada climate stations located in the Mackenzie Valley.

Climate station	MAAT (°C)	Annual precipitation (mm)	Annual snowfall (cm)
Tuktoyaktuk	−10.1	160.7	103.1
Inuvik	−8.2	240.6	158.6
Norman Wells	−5.1	294.4	161.5
Fort Simpson	−2.8	387.6	187.0

Note. MAAT indicates the use of the climate normal to calculate the mean annual air temperature. AMAT indicates the calculation of annual mean air temperature using less than 30 years of data (not the climate normal).

Figure 1. Site map with active layer monitoring sites. Sites are located both above and below treeline and monitor active layer thickness over both continuous and discontinuous permafrost. Sites with the longest records are located in the northern and southern ends of the delta. Permafrost zones from Heginbottom, Dubreuil, and Harker (1995) and ecoregion boundaries from Ecosystem Classification Group (2007, 2010, 2012).

Figure 2. Site photos with representative vegetation cover for each ecoregion sampled. More site details can be found in the Supplementary Material Table S1 (Smith et al. 2009; Duchesne, Riseborough, and Smith 2014).

Methods

Active layer monitoring sites throughout the study area were established to represent variable terrain conditions including surficial material and vegetation characteristics of the Mackenzie Valley (Nixon and Taylor 1994; Duchesne, Riseborough, and Smith 2014; Duchesne et al. 2015). Thaw tubes (e.g., Mackay 1973; S. L. Smith, Riseborough et al. 2009; S. L. Smith, Wolfe et al. 2009) were used to measure maximum annual thaw penetration and also the position of the ground surface at the time of maximum heave and ground surface subsidence. Maximum thaw depth is determined by the position of a glass bead that sits on the frozen surface and descends to the point of maximum thaw each year during the thaw season (Nixon and Taylor 1994). Measurements are made to 1 mm precision and are considered accurate to about 2 cm (Nixon and Taylor 1998). Thaw tubes are considered a superior technique compared to probing for determining active layer thickness because they do not require personnel to be present on the specific day of maximum thaw, which varies annually (Bonnaventure and Lamoureux 2013). Maximum thaw penetration and ground surface position at time of maximum thaw were used to calculate ALT. At some sites (n = 20) air and ground surface temperatures were recorded using single channel data loggers (Onset HOBO, accuracy ±0.4°C or Vemco, accuracy ±0.3°C) with air temperature sensors housed in radiation shields 1.5 m above the ground surface and ground surface temperature sensors buried at a depth of 3 to 7 cm (Duchesne, Riseborough, and Smith 2014).

Sites were classified by ecoregions by plotting the site location, using the specified coordinates on a location ecoregion map from Northwest Territories Environment and Natural Resources to assess classification. Of the fifty-five sites, thirty are in the tundra plains LAn ecoregion, eleven in the Taiga Plains HS ecoregion, six in the Taiga Plains LS ecoregion, four in the Taiga Plains MB ecoregion, and four in the Boreal and Taiga Cordillera (Table S1). Sites are not evenly distributed throughout the ecoregions partly due to logistical and accessibility issues.

Temporal gaps in the record were identified for each site. Most notable, there is a lack of recent data in the central portion of the valley, with most of the station records ending between 2008 and 2010. For each ecoregion, the annual minimum, maximum, and average of annual maximum ALTs were calculated. If there was only a record for one site in a given year for an ecoregion, no statistics were calculated, because it may have been misrepresentative of the ALT trends for that ecoregion.

Historical climate data and the Canadian monthly climate summaries for the nearest climate station (Environment Canada 2017) were used to determine monthly and yearly air temperature averages from 1990 to 2015 (the length of the monitoring record). At Tuktoyaktuk and Inuvik, multiple Environment Canada (EC) stations (n = 2 in close proximity, <400 m, with similar elevations, <0.3 m difference) were used to create a more continuous temperature record over the study period. A regression analysis was done for dates with data from both stations for comparison (Tuktoyaktuk r² = 0.67 and Inuvik r² = 0.99). Freezing and thawing degree-days were calculated for each year (September 1 to August 31) of the study period at each climate station through the summation of all the average daily temperatures below (FDD_a) and above (TDD_a) 0°C. FDD_a and TDD_a provide an indication of the duration and magnitude of freezing or thawing temperatures. At sites with on-site air and ground temperature measurements, FDD and TDD for the air (a) and ground surface (g) were calculated. Snowfall over the study period was compiled from the climate stations with data taken from several neighboring stations for Inuvik and Fort Good Hope. A regression analysis was done to compare the difference between the two stations (Inuvik r² = 0.57 and Fort Good Hope r² = 0.79). The regression analysis for Fort Good Hope was only used to infill missing snowfall data for April 2006. Missing snowfall values for the months of May to September were assumed to be 0 cm based on analysis of snowfall during this period over the entire 1990–2015 record and measurements of snowfall at other stations during the same month.

Air and ground surface FDD and TDD were used to calculate seasonal n-factors (ratio of surface temperature to air temperature; M. W. Smith and Riseborough 2002) for sites having on-site temperature data. N-factors are transfer functions representing the energy fluxes at a single location that can be used to estimate the surface offset (M. W. Smith and Riseborough 1996; Klene et al. 2001). Because the surface offset is influenced by vegetation and snow cover, which vary throughout the year, n-factors are typically divided into freezing (n_f) and thawing (n_t) n-factors, with n_f largely determined by snow cover in winter and the n_t characterizing vegetation cover in summer (M. W. Smith and Riseborough 2002). N_f and n_t were calculated using ground surface and air temperature data from October 1 of the previous year to April 30 of the current year and May 1 to September 30 of the current year, respectively (Equation (1)(1) ${\mathrm{n}}_f=\frac{FD{D}_g}{FD{D}_a}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{n}}_t=\frac{\mathrm{T}\mathrm{D}{\mathrm{D}}_g}{TD{D}_a}.$(1) ; M. W. Smith and Riseborough 2002; Duchesne, Riseborough, and Smith 2014).

(1) ${\mathrm{n}}_f=\frac{FD{D}_g}{FD{D}_a}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{n}}_t=\frac{\mathrm{T}\mathrm{D}{\mathrm{D}}_g}{TD{D}_a}.$(1)

Climatic data were examined using the downscaling program ClimateNA (Wang et al. 2016). This examined monthly reanalysis data (downscaled PRISM data) for each of the active layer monitoring sites were used to compare annual and seasonal air temperature trends over the study period. This allowed for the assessment of generalized climatic changes for each specific ecoregion to assess the relationship between air temperature and ALT changes. Additionally, a simple random forest and multivariable linear regression (MVLR) analyses were conducted to determine the relative importance of winter (FDD_a and FDD_g) and summer (TDD_a and TDD_g) conditions on ALT at sites with both ALT and temperature records. Random forest (R Studio 2020, randomForest) was run to determine variable importance ranking and MVLR (Excel, Analysis ToolPak) was completed to test each variable for statistical significance for all sites. Only MVLR was done for each ecoregion individually. These were selected for their straightforward and interpretable variable importance rankings (or statistical significance).

Results

Climate analysis

Records of air temperature from Environment Canada stations indicate a warming trend during the study period (1990–2015; Figure 3). At the northernmost station, Tuktoyaktuk (Figure 1), annual mean air temperature (AMAT) increased by 3.1°C from 1990 to 2015, whereas at the southernmost station, Fort Simpson AMAT increased by 2.7°C. This equates to warming trends of 0.12°C/decade and 0.10°C/decade at Tuktoyaktuk and Fort Simpson, respectively. The greatest warming (0.24°C/decade) occurred at Fort Good Hope. These values are all comparable to that estimated for the Mackenzie region in the most recent Arctic assessment (Arctic Monitoring and Assessment Programme 2017). Results of the ClimateNA reanalysis data indicate that warming is seasonally asymmetric with warming trends observed in autumn and winter (warming 0.8–1.3°C/decade and 1.0–2.0°C/decade with slope errors of 0.03 and 0.05, respectively) but no trend evident in spring and summer for all ecoregions (cooling 0.1–0.3°C/decade and warming 0–0.3°C/decade with slope errors of 0.05 and 0.03, respectively). Spatially, the average annual FDD_a decreased southward while the average annual TDD_a increased, with the Tuktoyaktuk station recording the highest FDD_a and lowest value of TDD_a and Fort Simpson station recording the lowest FDD_a and the highest TDD_a (Figure 4). Temporally, there was a slight decrease in FDD_a at each of the EC stations; however, the fit in each case was very weak (r² < 0.1 slope errors ranging from 8.7 to 28.5). The greatest temporal range in FDD_a occurred at the Wrigley station and the lowest occurred at the Tuktoyaktuk station (Figure 4). There was a slight but weak (r² < 0.1 slope errors ranging from 4.9–15.0) increase in TDD_a over the monitoring period at all EC stations. Over the monitoring period the greatest range in TDD_a was observed at the Wrigley station and the lowest occurred at the Inuvik station (Figure 5). Across every ecoregion with active layer monitoring sites, FDD_a varied more interannually than TDD_a (Figure 5). Spatial range in FDD_a within the ecoregions was greatest in the two northernmost ecoregions and smallest for the southern ecoregions (Figure 6). The same spatial pattern in variability was also observed for TDD_a, although the variation within each ecoregion was smaller than for FDD_a.

Figure 3. AMAT with the trendlines for Environment Canada climate stations (see Figure 1 for location) closest to active layer monitoring sites in the Mackenzie Valley (from Environment Canada 2017).

Figure 4. Average daily temperatures below (FDD_a) and above (TDD_a) for the EC stations in the Mackenzie Valley. Mean FDD_a and TDD_a for all years with data are represented by an X. For most stations, maximum and minimum FDD_a and TDD_a are shown by the ends of the whiskers. However, for some sites outliers are shown as circles after the whiskers and the ends of the whiskers instead show the value for one and a half times the interquartile range. The ends of the box show the first (25 percent) and third (75 percent) quartiles and the black line within the box shows the median. The years with the highest and lowest FDD_a and TDD_a are written respectively above and below the marker.

Figure 5. Average ALT for each ecoregion plotted with FDD_a/TDD_a and snowfall from the nearest EC station. (A) Tundra Plains LAn average ALT with FDD_a/TDD_a and snowfall from the Tuktoyaktuk EC station. (B) Taiga Plains HS average ALT with FDD_a/TDD_a and snowfall from the Inuvik EC station. (C) Taiga Plains LS average ALT with FDD_a/TDD_a and snowfall from the Norman Wells EC station. (D) Taiga Plains MB average ALT with FDD_a/TDD_a snowfall from the Fort Simpson EC station. (E) Boreal and Taiga Cordillera average ALT with FDD_a/TDD_a snowfall from the Fort Simpson EC station.

Figure 6. Box and whisker plot for FDD_a, TDD_a, FDD_g, and TDD_g measured at sites for each ecoregion. “g” is ground and “a” is air. X represents the mean and the circle indicates the value of outliers, which are points outside one and a half times the interquartile range.

Snowfall in the valley varies substantially annually and these fluctuations are not consistent between stations (Figure 7). For the study period (1990–2015), the lowest snowfall was recorded at Tuktoyaktuk (annual average 106 cm) with limited interannual variability (range 97 cm), whereas the most snowfall was recorded at Fort Simpson (annual average 199 cm), which also exhibited the largest interannual variability (range, 234 cm). Over the monitoring period trends in snowfall varied by station. At Tuktoyaktuk and Fort Simpson snowfall exhibited a positive temporal trend (r² = 0.37 and 0.31 with slope errors of 1.6 and 1.5, respectively). Snowfall seemed to decrease through time at Fort Good Hope (r² = 0.67 slope error of 1.3). Although there appeared to be trends in snowfall at Inuvik, Norman Wells, and Wrigley, the correlation was very weak (r² < 0.1 with slope errors of 2.5, 0.94, and 1.5, respectively).

Figure 7. Annual total snowfall data for Environment Canada stations in the Mackenzie Valley from 1990 to 2015 (Environment Canada 2017). Mean snowfall for all years with data is represented by an X. The maximum and minimum yearly snowfall at each station are represented by the whiskers. The black line within the box shows the median yearly snowfall at each station, and the ends of the box show the first (25 percent) and third (75 percent) quartiles. The years with maximum and minimum annual snowfall at each station are written above and below the marker, respectively.

Active layer monitoring

Regional

The average ALT for 1991 to 2014 ranged from less than 50 cm in the northern portion of the Mackenzie Valley in the Tundra Plains LAn ecoregion to greater than 160 cm in the southernmost part of the valley in the Taiga Plains MB ecoregion (Figure 8). Generally, the active layer is thinnest in the Tundra Plains LAn ecoregion and thickest in the Taiga Plains MB increasing moving south (Table 3 and Figure 9). However, this does not hold true for the Taiga Plains LS and Boreal and Taiga Cordillera ecoregions. Spatial range in ALT was greatest for the northern ecoregions and Taiga Plains MB. The spatial range in ALT for the Taiga Plains LS ecoregion is large due to one outlier site with a relatively thick active layer. The spatial range in ALT for the remaining sites is 42 cm.

Table 3. ALT metrics for each ecoregion.

Ecoregion	Average ALT (cm)	Spatial ALT range (cm)	Average temporal ALT range (cm)	Minimum ALT (cm)	Maximum ALT (cm)
TundraPlains LAn	67	145	22	26	171
TaigaPlains HS	92	132	27	49	181
Taiga Plain LS	86	128 (*60)	19	48	176
TaigaPlains MB	129	121	35	79	199
Boreal and	65	54	22	42	96
Taiga Cordillera

Note. Spatial range refers to the difference in the maximum and minimum ALT across all years and sites. Minimum and maximum refer to the thickest and thinnest active layers across all sites and years. Average temporal range refers to the average of the individual sites’ temporal range over the monitoring period. * indicates the spatial ALT for the ecoregion excluding the one outlying site.

Figure 8. Mean active layer thickness for all years and sites classified by ecoregion. Thirty monitoring sites are in the Tundra Plains LAn ecoregion, eleven are within the Taiga Plains HS ecoregion, six are in the Taiga Plains LS ecoregion, four are in the Taiga Plains MB ecoregion, and three are in the Boreal and Taiga Cordillera.

Figure 9. Mean and ranges of ALT by ecoregion along the Mackenzie Valley transect. Mean ALT for all years with data is represented by an X. The maximum and minimum ALT for all years and sites within an ecoregion are represented by the ends of the lines except for those with outliers (circles past the ends of the whiskers), in which case the ends of the whiskers represent the value one and a half times the interquartile range. The black line within the box shows the median annual ALT for each ecoregion, and the ends of the box show the first (25 percent) and third (75 percent) quartiles. The year and site with the thickest and thinnest active layers are written respectively above and below the marker. Some sites maintained the thickest or thinnest active layer for multiple years, and for these sites the years indicate this duration. More details on the individual stations can be found in the Supplementary Material in Table S1.

Temporally, sites with thinnest active layers generally remained thin throughout the monitoring period. Similarly, sites with the thickest active layers remained the thickest. Over the monitoring period ALTs in the Taiga Plains HS and MB ecoregions varied the most through time while ALTs in the Tundra Plains LAn, Taiga Plains LS, and Boreal and Taiga Cordillera showed the least variability.

Air–ground surface connectivity also varied spatially and temporally between ecoregions. Average n_f values for each ecoregion decreased from north to south, with the highest mean n_f (0.43) in the Tundra Plains LAn ecoregion and the lowest value (0.27) in both the Taiga Plains MB and Boreal and Taiga Cordillera ecoregions (Table 3 and Figure 10). The two northernmost ecoregions (Tundra Plains LAn and Taiga Plains HS) had the largest spatial range in n_f, whereas the Taiga Plains MB and Boreal and Taiga Cordillera had the smallest range. Temporally, variability in n_f over the monitoring period was similar for each ecoregion (0.2–0.12) but was lowest for the two northern ecoregions and highest for the Taiga Plains MB and Boreal and Taiga Cordillera ecoregions. Within each ecoregion n_f showed more temporal variation at sites with higher n_f values. Average n_t was more similar between ecoregions than n_f. Like n_f, the spatial range in n_t was highest for the northern ecoregions (Tundra Plains LAn and Taiga Plains HS) and smallest for the Taiga Plains MB and Boreal and Taiga Cordillera ecoregions (Table 4 and Figure 10). The range in n_t interannually was greatest for the Tundra Plains LAn and Boreal and Taiga Cordillera and smallest for the Taiga Plains MB ecoregion.

Table 4. Average n_f and n_t values and ALT for sites with data.

Site	Ecoregion	Average nf	Average nt	Average ALT (cm)	FDDg (°C days)	TDDg (°C days)
90-TT-04	Tundra Plains LAn	0.45	0.33	87	2,033	428
90-TT-05	Tundra Plains LAn	0.57	0.53	75	4,411	524
90-TT-06	Tundra Plains LAn	0.37	0.28	48	1,697	325
90-TT-13	Tundra Plains LAn	0.64	0.58	63	2,616	386
91-TT-C	Tundra Plains LAn	0.09	0.98	125	403	801
92-TT-01	Tundra Plains LAn	0.49	0.29	47	2,215	311
92-TT-02	Tundra Plains LAn	0.37	0.47	56	1,693	415
90-TT-17	Taiga Plains HS	0.09	0.71	147	421	997
91-TT-12	Taiga Plains HS	0.55	0.70	73	2,288	771
91-TT-13	Taiga Plains HS	0.24	0.89	134	997	1,085
91-TT-14	Taiga Plains HS	0.29	0.59	109	1,277	782
91-TT-15	Taiga Plains HS	0.32	0.57	77	1,482	890
91-TT-16	Taiga Plains HS	0.30	0.57	89	1,332	863
91-TT-22	Taiga Plains LS	0.25	0.41	62	992	688
92-TT-04*	Taiga Plains MB	0.29	0.61	124	1,001	1,202
92-TT-05	Taiga Plains MB	0.22	0.55	184	816	1,102
92-TT-06	Taiga Plains MB	0.22	0.50	137	654	978
92-TT-07	Taiga Plains MB	0.34	0.73	86	1,191	1,455
92-TT-08	B&T Cordillera	0.24	0.49	79	916	913
92-TT-10	B&T Cordillera	0.29	0.50	64	995	961

Note. *For site 92-TT-04 the air temperature data from site 92-TT-05 was used due to the proximity of the two stations.

Figure 10. Box and whisker plot for the freezing n-factor (n_f) and the thawing n-factor (n_t) for each ecoregion. X represents the mean and the circle indicates the value of outliers (one and a half times the interquartile range).

Spatial and temporal variability in FDD_g and TDD_g was also observed between ecoregions. FDD_g generally decreased from north to south (Figure 6). Northern ecoregions also had the largest spatial variability in FDD_g between sites. Interannually, the greatest range in FDD_g occurred in the Taiga Plains HS ecoregion and the lowest was in the Taiga Plains MB ecoregion. Average TDD_g exhibited a positive spatial trend from the northern ecoregions to the southern except for the Taiga Plains LS ecoregion, which only had one site. Temporally, the Boreal and Taiga Cordillera ecoregion showed the greatest variability in TDD_g, whereas the Taiga Plains MB showed the least over the study period. Temporal trends for n_f, n_t, FDD_g, and TDD_g were difficult to determine due to data gaps and limited years with data.

Variable importance ranking in the random forest analysis identified FDD_g, TDD_g, and n_f as the three most important factors determining ALT (in this order), followed by n_t, FDD_a, and, lastly, TDD_a. For overall average ALT the MVLR analysis produced an r² of 0.61, indicating that only FDD_g was significant (p = .032). Using all measurements across all sites and years, FDD_g, TDD_g, and TDD_a were all significant (p < .001 for FDD_g and TDD_g and p = .006 for TDD_a) with an r² of 0.45. When divided by ecoregion, the MVLR analysis showed variability in significant variables between ecoregions. For the Tundra Plains LAn ecoregion only TDD_g was significant (p = .001, r² = 0.55), whereas for the Taiga Plains MB ecoregion only FDD_g was significant (p = .03, r² = 0.46). The Taiga Plains HS and Boreal and Taiga Cordillera MVLR analysis showed FDD_g and TDD_g to be significant with p < .001 and p = .001 (r² = 0.83) for Taiga Plains HS and p = .001 and 0.047 (r² = 0.46) for the Boreal and Taiga Cordillera. The results indicate the importance of considering antecedent winter conditions when analyzing ALT, especially in the southern portion of the valley where they are shown to have a higher level of significance.

Site-specific

Site-specific changes in ALT and measured temperature parameters showed unique temporal changes over the duration of the study. At sites 91-TT-13 (Taiga Plains HS), 91-TT-22 (Taiga Plains LS), 92-TT-04 (Taiga Plains MB), and 92-TT-10 (B&T Cordillera), FDD_g and n_f showed a degree of variability through time (Figure 11). Of these sites, 91-TT-13, 92-TT-04, and 92-TT-10 also showed increasing ALT through time (0.5, 1.7, and 0.5 cm yr⁻¹, respectively). At site 92-TT-04 the increase in ALT over the study period corresponded with a drop in FDD_g and n_f as well as an apparent increase in n_t and a large degree of interannual variation in TDD_g. Sites 91-TT-22 and 92-TT-10 had limited ALT response and showed contrasting warming winter and cooling summer ground surface temperature trends.

Figure 11. Average ALT, TDD_g, and n_t and FDD_g and n_f for the preceding year. Slope and slope error for each series are given. (A) Site 90-TT-04 in the Tundra Plains LAn ecoregion. (B) Site 91-TT-13 in the Taiga Plains HS ecoregion. (C) Site 91-TT-22 in the Taiga Plains LS ecoregion. (D) Site 92-TT-04 in the Taiga Plains MB ecoregion. (E) Site 92-TT-10 Boreal and Taiga Cordillera.

Discussion

It is important to note that the spatial variability was greater within ecoregions than between ecoregions likely due to the variations in snow cover, vegetation, and substrate within each ecoregion that control ALT on a subecoregional scale (Bonnaventure and Lamoureux 2013). This highlights the difficulty in using general bioclimatic variables to broadly characterize ALT, particularly in areas dominated by shrubs. Given the high degree of spatial variability in shrub-dominated environments, an ecoregional approach may be inappropriate especially at a local scale. However, for broad-scale modeling studies of ALT, the ecoregional approach may still be advantageous because unique spatial and temporal patterns can still be discerned. This is largely a product of variation in the dominant influence on ALT between each ecoregion and will be explored further in the following sections.

FDD and TDD

The random forest and MVLR highlighted the importance of both winter and summer conditions on ALT. In previous studies the impact of winter conditions has been largely overlooked with an emphasis placed on using summer conditions, primarily TDD_a, to infer about the state of the active layer (Shiklomanov and Nelson 2002; Mazhitova et al. 2004). However, in recent years the importance of antecedent winter conditions has been documented (Wilcox et al. 2019; Grünberg et al. 2020). This analysis reaffirms the importance of considering antecedent winter ground conditions in regional ALT studies, especially in more southern environments where those were shown to be most important.

Moving north to south, average ALT for ecoregions generally increases corresponding to a decrease in FDD_a and increase in TDD_a (Figures 5 and 9). The exception to this trend is the average ALT in the Boreal and Taiga Cordillera (68 cm), which is similar to that of the Tundra Plains LAn (67 cm) despite the difference in FDD_a (2954°C days in the Cordillera compared to 4480°C days in the Tundra Plains) and TDD_a (2132°C days in the Cordillera compared to 934°C days in the Tundra Plains). However, the ground surface thermal conditions for each ecoregion show some similarity during the freezing (FDD_g of 955°C days for the Cordillera and 1880°C days in the Tundra Plains) and thawing seasons (TDD_g of 937°C days in the Cordillera and 456°C days for the Tundra Plains). The similarity in ALTs and the ground thermal regime (particularly the summer conditions) despite large differences in the air thermal regime is likely indicative of in influence of protective vegetation cover or organic material (Shur and Jorgenson 2007). Despite similar air thermal conditions, the difference in the average ALT (58 cm) for the two southernmost ecoregions (Boreal and Taiga Cordillera and Taiga Plains MB) was greater than the difference between the Boreal and Taiga Cordillera and the Tundra Plain LAn ecoregions (1 cm). This difference between the ALT of the two southernmost ecoregions is likely the result of mosses and organic material at sites in the Boreal and Taiga Cordillera that allow for the thin active layer despite less favorable air temperature (Bonnaventure and Lamoureux 2013). However, there are only a small number of sites with data in the Boreal and Taiga ecoregion, which may not be representative of the entire ecoregion.

Vegetation and snow cover dynamics

The impact of snow cover on the active layer and permafrost on a yearly scale is dependent on the timing, duration, thickness, structure, and physical and thermal properties of the snowpack (Hinkel and Nelson 2003; Ling and Zhang 2003; Zhang et al. 2005). The accumulation of snow cover in autumn and early winter is shown to have a strong influence on ground temperature by insulating the ground from the low winter air temperatures, resulting in warmer ground surface temperatures compared to locations with thin or no snow cover (Gold 1963; Goodrich 1982; Zhang 2005). N_f and FDD_g values are often used to show air–ground surface winter connectivity, which is largely influenced by snow conditions at a site. Intersite variability in vegetation characteristics is one of the likely reasons for variability in ALT within the ecoregions and is consistent with other local studies (Fisher et al. 2016; S. L. Smith et al. 2016). Based on the variability in n_t data this appears to partly explain a reason for the spatial variability between and within ecoregions.

A previous ecoregional assessment of n_f in the Mackenzie Valley found that n_f generally increased moving northward, largely indicating an increased connectivity between the air and ground surface at northern sites representing less and more dense snow cover (S. L. Smith et al. 2016). When comparing ecoregions, the same trend is apparent with the northernmost Tundra Plains LAn ecoregion having the highest mean n_f and the southernmost Taiga Plains MB and Boreal Cordillera ecoregions having the lowest mean n_f (Figure 10). Within ecoregions lower n_f values are associated with greater ALT because these sites experience warmer winter ground surface conditions than sites with high n_f. Additionally, the variability in n_f was highest in the Tundra Plains LAn ecoregion (Figure 6), which was attributed to highly variable spatial snow cover within the ecoregion, likely due to the presence of shrubs at some sites which trap and retain blowing snow (Yang and Woo 1999; Sturm et al. 2001; Zhang 2005; Palmer, Burn, and Kokelj 2012; Bonnaventure, Lamoureux, and Favaro 2017; Garibaldi, Bonnaventure, and Lamoureux 2021). This combination of vegetation and snow cover may also explain why shrub-dominated ecoregions (Tundra Plains LAn and Taiga Plains HS) showed the greatest variability in ALT within ecoregions. Forest-dominated ecoregions likely show less interecoregional variability due to more uniform snow conditions. The greatest observed ALT and n_f and n_t variability in the Tundra Plains LAn ecoregion may also just be a result of a higher level of sampling sites (n = 30) compared to any other ecoregion. This allowed for a wider range of site characteristics to be sampled, giving this ecoregion the most robust variability of sites. Changes in site-specific ALT may also reflect changes in vegetation and snow cover dynamics.

Summer temperatures, represented by TDD_a, TDD_g, and the corresponding surface offset, n_t, were also shown to influence ALT. Between ecoregions, n_t decreased from the northern ecoregions to the southern ecoregions in response to changing vegetation conditions. n_t is highest for the Tundra Plains LAn ecoregions due to the lack of shading during summer as a result of no tree cover and limited vegetation or topographic shadowing. This results in less uniform heating between the ground and the air and can be amplified when darker surface soils are present. In the southern forested ecoregions, n_t is low due to shading effects of vegetation and also thicker organic layers muting the response of the ground temperatures to air temperatures through increased evaporation and differential thermal conductivity (Shur and Jorgenson 2007; Fisher et al. 2016).

At the northernmost sites, especially those in the Tundra Plains LAn ecoregion including 90-TT-04, there appears to be more of a direct link between FDD_g and ALT as changes in some years (2007–2009) closely follow each other (Figure 11). Additionally, there was no apparent negative trend in ALT despite a decrease in TDD_g (Figure 11). The influence of increases in TDD_g was amplified moving southward as ALT responded more to changes in TDD_g alone and in conjunction with changes in FDD_g such as at sites 91-TT-22 and 92-TT-04 (Figure 11). Additionally, the warming of the ground surface in winter at sites (decrease in FDD_g) in the Taiga Plains MB ecoregion may reflect the increase in the overall snowfall of the area because the Fort Simpson climate station indicates about a 7 cm increase in snowfall per year over the study period (Chapin et al. 2005).

Because ALT is influenced by snow accumulation and vegetation, a disturbance or change in vegetation may result in changes to ALT through changes to snow entrapment or loss of shading, transpiration, and organic matter (Shur and Jorgenson 2007; Fisher et al. 2016; Way and Lapalme 2021). This was observed at site 92-TT-12 in the Taiga Plains LS ecoregion, which burned in a fire in 1995. In the twelve years following the fire, ALT increased by approximately 44 cm. This magnitude of change during this time frame is consistent with other studies examining active layer changes postfire (S. L. Smith, Riseborough, and Bonnaventure 2015; Holloway and Lewkowicz 2020; Holloway et al. 2020). At this site, the active layer thickening is likely the result of surface modification during the fire because no other site nearby or in the ecoregion shows this trend. An additional contributing factor to active layer development could be the presence of surface water that does not freeze to the ground surface, which can substantially keep the ground warm during winter, resulting in an overall warming effect (Harris 2002; Smith, Wolfe et al. 2009). This is likely the case at site 91-TT-C, which has been observed to occasionally flood.

Substrate

Substrate, particularly grain size, moisture, and ice content moderate ALT through latent heat exchanges and the differential thermal conductivity between the frozen and thawed state (S. L. Smith, Burgess, and Nixon 2001; Shur and Jorgenson 2007; O’Donnell et al. 2009; Bonnaventure and Lamoureux 2013). Surface subsidence results as thaw progresses deeper into ice-rich ground. Seasonal freeze–thaw of the active layer results in latent heat exchanges and seasonal heave and settlement. However, ongoing settlement can occur as the thaw progresses deeper over time, thawing new ice-rich material (O’Neill, Smith, and Duchesne 2019). Therefore, at sites with high excess ice content, ALT can remain relatively stable as melting of ground ice results in surface subsidence (Duchesne et al. 2015; O’Neill, Smith, and Duchesne 2019). Sites with thin active layers and very little annual variability are often associated with ice-rich material susceptible to thaw settlement.

The results indicate that sites with frost-susceptible substrate such as ice-rich silts and clay are often associated with higher subsidence than variability in ALT (Figure 12). Conversely, sites with ice-poor sandy substrates often showed little subsidence compared to variability in ALT over the study period (Shur, Hinkel, and Nelson 2005). There was no pattern relating subsidence and ecoregion because the type of substrate and ice content varies widely within ecoregions, likely contributing to the intra-ecoregion ALT variability.

Figure 12. Comparison of ALT range and subsidence for each site classified based on substrate texture (Geological Survey of Canada data files updated from S. L. Smith, Riseborough et al. 2009). Subsidence was measured cumulatively over all years of observations. Range of ALT refers to the difference between the maximum and minimum ALT for all years at a site. Red symbols indicate the average subsidence for that substrate.

Some of the sites can be divided into four classes based on the range in ALT and subsidence (low ALT range high subsidence, low ALT range low subsidence, high ALT range low subsidence, high ALT range high subsidence). Examples of low ALT variability and high subsidence sites include 91-TT-11, which has a thin active layer (66 cm) with very little annual variability but with a surface subsidence of 25 cm (over nineteen years), and site 91-TT-02, which has similar active layer conditions and surface subsidence of 36 cm over seventeen years (Geological Survey of Canada data files updated from Smith, Riseborough et al. 2009). These sites both have fine-grained clay substrate, are low elevation, are flat, and are located in an area that has a medium potential for segregated ground ice content according to recent ground ice maps (O’Neill, Wolfe, and Duchesne 2019). At these sites, the high subsidence and consolidation of the thaw layer result in smaller changes in ALT. Sites with low ALT variation and low subsidence include 90-TT-01 and 90-TT-12; both have fine-grained, silt substrate with high ground ice content limiting the progression of thaw (O’Neill, Wolfe, and Duchesne 2019). These sites are two of the northernmost sites with two of the thinnest active layers. Therefore, it is possible that the low subsidence at these sites is a result of limited thaw penetration into the ice-rich layers of the substrate. A high range in ALT with low subsidence occurs at sites 90-TT-06, 91-TT-13, and 92-TT-05 likely due to limited ground ice. Site 90-TT-06 shows 6 cm subsidence over twenty years with a 30 cm range in ALT. Though this site has a small average ALT (48 cm), it has a sand substrate in which excess ice is unlikely to be present. Site 91-TT-13 is ice-poor silt, which may explain the limited subsidence (8 cm), large variation in ALT (39 cm over twenty-two years), and thick average active layer (134 cm). Lastly, sites with both high ALT variation and high subsidence such as 92-TT-03 and 92-TT-12 are likely in places where ice-rich material may have already thawed and consolidated early in the record. At 92-TT-03 little change in ALT (11 cm) and substantial subsidence (26 cm) occurred from 1994 to 2004, followed by larger changes in ALT (34 cm) and minimal subsidence from 2004 to 2014. This site is classified as ice-rich silt. A similar pattern occurs at 92-TT-12 with limited variation in ALT (28 cm) and substantial subsidence (43 cm) during the initial part of the record (1993–2001) followed by higher variability in ALT (up to 10 cm fluctuations between years) and minimal subsidence (8 cm). This site is classified as ice-rich peat. Higher ice contents are commonly found at the top of the permafrost (Shur, Hinkel, and Nelson 2005), and the larger subsidence observed earlier in the record is likely the result of thawing of this ice-rich material. As thaw progresses into material containing less ice, less subsidence occurs and larger variation in ALT is observed. Latent heat effects also decrease as thaw progresses into the lower ice-content material (Shur, Hinkel, and Nelson 2005).

Temporal trends in ALT

Potentially a key rationale for utilizing an ecoregional approach is the ability to assess broad temporal trends in ALT based on bioclimatic variables hypothesized to impact the response of permafrost to climate change (Shur and Jorgenson 2007). Using the three-year running means for all annual ALT in each ecoregion over the study period, only the Taiga Plains HS showed a strong positive trend in average ALT (r² = 0.8) despite the positive air temperature trend in all regions (Figure 13). ALT for the Tundra Plains LAn and Boreal and Taiga Cordillera showed a slight positive trend with low r² values of 0.1 and 0.4, respectively. The Taiga Plains LS and Taiga Plains MB ecoregions showed a negative trend or no trend with low r² values of 0.2 and 0.0, respectively. The trends may be impacted by the number of sites per ecoregion because an ecoregion with a smaller number of sites is more influenced by changes in individual sites than an ecoregion with a larger number of sites. In this case the trend or lack thereof would not be representative of the entire ecoregion. This is more likely to be the case in the Taiga Plains LS and MB ecoregions because they have fewer sites than the Taiga Plains HS and Tundra Plains LAn ecoregions and have one site with a substantially different ALT than the other sites in the ecoregion. The variation in ALT responses to similar warming throughout the valley may demonstrate the variable sensitivity to climatological or ecological disturbances depending on the ecoregion. For example, permafrost in most tundra environments is considered to be climate driven (e.g., Shur and Jorgenson 2007), and ALT in these ecoregions is a product of low air temperatures and is thus particularly susceptible to long-term changes in climate. This is due to a limited ecological buffer shielding the ground from the direct impacts of increasing air temperatures. The thickening of the active layer in the Taiga Plains HS ecoregion may demonstrate this however, the Tundra Plains LAn ecoregion did not show this same thickening, possibly as a result of higher ground ice (O’Neill et al. 2019). Other forested ecoregions with well-developed organic surface layers including those in the Boreal Cordillera can be classified as climate-driven, ecosystem-modified, and ecosystem-protected permafrost. Thus, the active layer is potentially less sensitive to short- or medium-term changes in air temperature but is more sensitive to vegetation disturbance including fire and hydrological and anthropogenic changes (Shur and Jorgenson 2007; Fisher et al. 2016). Therefore, it is unlikely that even on the time frame of this study or into the short-term future these locations will display a substantial response to climatic change through active layer thickening unless paired with substantial vegetation disturbance or pond or wetland expansion (Holloway and Lewkowicz 2020). This may explain why the ecoregions with more developed vegetation cover did not on average show a response over the monitoring period. As a result, the variability in response of the active layer depends on the type of disturbance and thus it is important to understand this in order to properly conceptualize the timing of associated landscape hazards, which may result from permafrost degradation.

Figure 13. Three-year running means for average active layer thickness for each ecoregion. R², slope, and slope error are given for each series.

Summary and conclusions

ALT at monitoring sites throughout the Mackenzie Valley exhibits a great deal of spatial and temporal variability within and across ecoregions because ALT is a product of complex interactions between climate, vegetation, snow cover, and soil type. The Tundra Plains LAn ecoregion showed the largest spatial range in ALT, whereas the Taiga Plains MB ecoregion has the greatest average ALT. Thick active layers at certain sites in the Tundra Plains LAn and Taiga Plains HS ecoregions likely result from the increased accumulation of snow by willow shrubs, resulting in warmer ground conditions than tundra upland sites that are wind scoured and have thin active layers. This causes large ranges in ALT in these ecoregions. Generally, air temperature, including winter conditions, is the major influence because average ALT increases moving southward, a result of air temperature increase and greater ALT occurring in years with warmer winters and summers. However, the relative importance of winter and summer thermal conditions on ALT depends on the ecoregion. The characteristic vegetation type of an ecoregion may also be used as an explanation for the average ALT; however, this is more difficult to use because dominant vegetation may vary between sites within an ecoregion. Overall, winter conditions can have significant influence on ALT.

ALT at some sites was less responsive to shorter-term variation in air temperature, potentially due to the moderating effect of vegetation, organic layer, and snow. For these locations, disturbance, rather than changing air temperature, may have a larger impact on ALT due to the surface modification. Where ice-rich permafrost is present, increases in summer thaw due to warming and increased disturbance will likely result in increased ground subsidence. This subsidence may result in landscape changes that impact hydrology and drainage as well as ecosystems throughout the valley, potentially causing further changes in ALT, and increase permafrost degradation with increased ponding. Understanding the relation between ALT, climate, dominant vegetation type, and ecoregion is essential to assess how ALT will respond to changing air temperatures, snowfall, and vegetation under climate warming.

Acknowledgments

Logistical support was provided by Polar Continental Shelf Program and Inuvik Research Institute. We are thankful for the support and assistance provided by communities in the Northwest Territories. Numerous colleagues have contributed to data collection over the years. We also thank A. Leblanc for review prior to submission. Data management and processing support was provided by J. Chartrand.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental material for this article can be accessed on the publisher’s website.

Additional information

Funding

This work was supported by the Government of Canada, Natural Sciences and Engineering Research Council of Canada.