Highlights

  • Rapid shrubification has occurred in the Torngat Mountains, Nunatsiavut, Canada

  • Field observations indicate the area is characterized by widespread but discontinuous permafrost

  • Local differences in plant size and cover are moderated by the soil moisture regime

Introduction

Rapid changes have occurred in the distribution (Tape and others 2006; Myers-Smith and Hik 2017) and abundance (Sturm and others 2001b; Elmendorf and others 2012b) of Arctic shrubs in recent decades as climatic limitations on their establishment, growth, and productivity have eased. Large areas of tundra have recently undergone a greening trend (Ju and Masek 2016); however, responses to climate change have been heterogeneous across space and time (Elmendorf and others 2012b; Myers-Smith and others 2020) due to the interconnectedness of Arctic vegetation with various biological and environmental feedbacks (Chapin III and others 2005; Sturm 2005; Morrissette-Boileau and others 2018). Arctic plants provide essential resources and habitat for animals, influence biogeochemical cycling, and impact how humans interact with and make use of the land (Post and others 2009; Pearson and others 2013). Thus, understanding how changes in vegetation characteristics manifest across spatial and temporal scales is critical for predicting future impacts of climate change on tundra ecosystems and the people who use them.

Of the plants that grow beyond the forest-tundra ecotone, shrubs disproportionately influence local environments because of their relatively large physical presence compared to other, low-statured lifeforms. Changes in the size, distribution, and density of shrub species, termed ‘shrubification’, have been widely observed throughout the Arctic using repeat photography (Tape and others 2006; Lantz and others 2013; Ropars and others 2015b), field experiments (Ackerman and others 2017; Myers-Smith and Hik 2017), and satellite remote sensing (Fraser and others 2011; McManus and others 2012). Shrubs are able to respond quickly to improved conditions through vegetative reproduction and growth (Chapin III and others 1995; Zamin and Grogan 2012), and are often associated with declines in moss and lichen species through shading and litter production (Chapin III and others 1995; Cornelissen and others 2001; Joly and others 2009; Chagnon and Boudreau 2019). Such changes in plant community dynamics have cascading implications for local animal populations; for instance, shrubification may provide new habitat for bird species (Mizel and others 2016; Whitaker 2017), but may be detrimental for caribou which primarily forage for lichen in the winter (Joly and others 2009; Schmelzer and others 2020). Caribou in particular can have a dynamic influence on shrub establishment and growth (Andruko and others 2020) as browsing pressure is linked to overall herd size.

Changes in shrub cover and height also indirectly influence tundra landscapes through various feedback mechanisms, the most impactful of which may be the connection between shrubs and snow redistribution. Upright shrubs trap blowing snow in their dense network of branches, altering local snow accumulation (Pomeroy and others 1997; Sturm and others 2001a) and insulating the ground throughout the winter (Myers-Smith and Hik 2013; Roy-Léveillée and others 2014; Domine and others 2016; Paradis and others 2016). It is hypothesized that winter insulated soils could further amplify shrub growth through elevated rates of organic matter decomposition (Sturm and others 2005), as plants in tundra environments are often nutrient-limited (Weintraub and Schimel 2005). Changes to soil temperatures also affect and are affected by, near-surface permafrost. At local scales, thawing permafrost may enhance shrub establishment by exposing mineral soil (Lantz and others 2009; Frost and others 2013) and altering soil moisture regimes (Osterkamp and others 2009). At regional scales, increasing quantities of terrestrial carbon from peatland and tundra soils are expected to be released due to permafrost thaw in the coming decades (MacDougall and Knutti 2016; Turetsky and others 2020). It has been suggested that the enhanced shading of taller shrubs during the summer could have a cooling influence on permafrost (Blok and others 2010); however, studies have shown the magnitude of this effect varies (for example, Johansson and others 2013; Myers-Smith and Hik 2013; Paradis and others 2016) and is generally overwhelmed by the larger magnitude of winter warming through enhanced snow cover (Zhang 2005; Lawrence and Swenson 2011; Loranty and others 2011, 2018; Domine and others 2016).

It is clear that multiple biotic and abiotic processes and feedback loops drive vegetation change in northern environments (Myers-Smith and others 2011,2020). However, much of our knowledge is based on research from western North American Arctic and alpine tundra ecosystems. The eastern Arctic in North America has a stronger maritime climate, was more recently deglaciated (Dyke 2004), began warming later than north-western Canada (Turner and others 2007), and is strongly influenced by the effects of sea ice cover (Banfield and Jacobs 1998). These climatic differences raise the possibility that the direction, rate, or strength of interactions between plants and their environment may also differ. Remote-sensing data show that the eastern Canadian Arctic and Subarctic have recently undergone some of the highest rates of vegetation greening in North America (Fraser and others 2011; Ju and Masek 2016), but there are few in situ studies demonstrating this trend, with most of these studies from lowland sites in Nunavik and eastern Québec (for example, Tremblay and others 2012; Provencher-Nolet and others 2014; Ropars and others 2015a; Paradis and others 2016; Morrissette-Boileau and others 2018). The Torngat Mountains in northern Labrador comprise the southern limit of both the Canadian Arctic and the Arctic Cordillera (Ponomarenko and McLennan 2010; Riley and others 2013), and represent a critical threshold between Subarctic and Arctic ecosystems. Recent climate change in northern Labrador has affected wildlife (Whitaker 2017), glaciers (Way and others 2015; Barrand and others 2017), and terrestrial ecosystems (Siegwart Collier 2020). In the northern Nunatsiavut region of Labrador, sea ice has declined by 33% per decade from 1981 to 2010 (Barrette and others 2020), locally amplifying the influence of global climate warming. Inuit in Nunatsiavut and Nunavik have observed notable changes in the biophysical environment, including weather, snow, and ice conditions (Downing and Cuerrier 2011; Rapinski and others 2018) and in plant communities (Parks Canada 2008; Cuerrier and others 2015; Siegwart Collier 2020). These many environmental changes have impacts on the health, well-being, and sense of place of Inuit in Labrador (Cunsolo Willox and others 2012). The Torngat Mountains are expected to undergo considerable ecosystem modification over the next century due to changes in the region’s climate (Barrette and others 2020) and proximity to the forest-tundra ecotone. An integrated understanding of relationships between vegetation change and the broader environment in the Torngat Mountains remains elusive due to difficulty in accessing this remote and isolated mountain range, and the interdisciplinary methods required to study regional environmental change.

Here, we present the results of ten years of multidisciplinary research on vegetation–environment interactions within the Nakvak Brook watershed in Torngat Mountains National Park of Canada (TMNP; Tongait KakKasuangita SilakKijapvinga in Inuktitut). The area is an Inuit homeland that falls within an overlapping area between the Labrador (that is, Nunatsiavut) and Nunavik Inuit Land Claims Agreements, and the park is co-operatively managed with Inuit from Nunatsiavut and Nunavik. The study was motivated by questions regarding ecosystem change raised during more than a decade of discussions with the Torngat Mountains National Park Cooperative Management Board, as well as time spent sharing observations on the land with Inuit [including Park staff, members of the Cooperative Management Board, local Inuit students, professional bear guards, and a co-author who is of Inuit descent (R.G. Way)]. We combine long-term remote-sensing data with field-based observations of vegetation change and permafrost to characterize the climatic and biophysical environment, and to assess the sensitivity of ecological and geomorphological systems to ongoing climate change. Combining field observations with time series data allows us to contextualize local conditions within a recent history of regional environmental change in a remote, topographically complex, marine-influenced mountain range, and offers insight into how further changes might manifest in the future.

Study Location

The study area at Nakvak Brook (Pitukkik in Inuktitut) is situated in a coastal Cordilleran landscape at the transition between the Canadian low Arctic and Subarctic regions in northern Labrador (58.64° N, 63.35° W; 420 m above sea level (asl); Figure 1). The nearby Nakvak Brook valley leads to the Koroc River valley, Québec, and has been used by Inuit as a travel route between the Labrador coast and Ungava Bay for generations. Atmospheric reanalysis data (ERA5; Hersbach and others 2020), statistically downscaled to the location and elevation of the site, indicate that mean temperatures of the coldest and warmest months are − 25 °C (February) and 8 °C (August), respectively, and that the area experiences a mean annual precipitation of 875 mm, much of which accumulates as snow from October to May (1989–2018 climate normals). The regional climate is influenced by the Labrador Current, which promotes cold maritime conditions during the summer growing season and delays peak summer warmth to August, and by seasonal sea ice cover, which causes the climate in the region to become increasingly continental during winter (Banfield and Jacobs 1998).

Figure 1
figure 1

A Map of the study region in Torngat Mountains National Park (TMNP), northern Labrador. The red square indicates the study area at Nakvak Brook. B Spatial extent of the NDVI imagery used in this study; the red box indicates the extent of plot-scale measurements. C Schematic of the study site showing locations of vegetation plots, three permafrost electrical resistivity tomography/shrub transects, and water features.

Geomorphologically, the Nakvak Brook research site is a relatively flat upland basin surrounding a shallow lake, with relief to the east that rises to 900 m asl. The study site is located about 25 km inland from the coast of the Labrador Sea and ca. 15 km from the outflow of Nakvak Brook in Saglek Fiord. Local high points have relatively minimal (inferred as < 1 m thick) surficial cover, whereas riparian deposits are thicker (> 3 m) and range from coarse to medium in particle size. The finest soils exist at the lowest elevations next to streams entering and leaving the largest central lake (Figure 1C). Soil moisture is highly variable within the study site with a network of small surface runoff channels and alternating areas of saturated and dry ground. Visual observations (L. Hermanutz) suggest that wet areas have undergone a drying trend in the last decade, supported by limited soil moisture measurements at the site (Supplementary Figure 1). According to the Permafrost Map of Canada (Heginbottom and others 1995), the region is underlain by extensive discontinuous to continuous permafrost; however, recent studies classify the broader region as extensive discontinuous permafrost (Way and Lewkowicz 2016; Obu and others 2019).

Ground cover at Nakvak Brook is primarily composed of low tundra vegetation and exposed rock. Arctic dwarf birch (Betula glandulosa Michx. (AvâlaKiak in Inuktitut)), willow (Salix herbacea L.; S. arctica Pall (UKaujak); Salix spp.), and blueberry (Vaccinium uliginosum L. (Kigutanginnak)) are the dominant deciduous shrub species in the area. Dry areas tend to have a greater abundance of lichen and evergreen shrubs, whereas mosses and graminoids are more common in wet areas. The larger TMNP region, including the Nakvak Brook watershed, has experienced a dramatic increase in vegetation and decrease in bare ground over the past several decades, including rapid expansion of tall shrubs in valley bottoms (Fraser and others 2011; Quirouette and Zorn 2015; Ju and Masek 2016). At the same time, there has been a significant decline in number of caribou in the Torngat Mountains since the 1990s, in large part due to the collapse of the migratory George River population, which used to range into the region, and possibly also due to a decline in the locally resident, but much smaller, Torngat Mountains population (Couturier and others 2015; COSEWIC 2017). The size of the Torngat Mountains population has recently increased slightly (ca. to 1300 individuals in 2017), but numbers remain low and the range of the remnant George River population no longer includes this region (Couturier and others 2018). Thus, it can be presumed that caribou browsing pressure has declined during the past 20 years.

Methods

Low Altitude Imagery From Remotely Piloted Aircraft

Low altitude remotely piloted aircraft (RPA) (DJI Phantom 4 Pro quadcopter) imagery (N = 1068) was acquired at Nakvak Brook on 4 August 2018 and was used to characterize the local geomorphological context and spatial distribution of surface features (large vegetation patches, streams, lakes, frost sorted features) in relation to field data. RPA imagery and a description of collection and processing methods are provided in Supplementary Figure 2.

Decadal-Scale Changes in Air Temperature and Vegetation

Climate and surface reflectance data were acquired to generate local records of temperature and vegetation change for the study area. Prior research in Labrador has shown that scaling atmospheric reanalysis and gridded climate products can improve local temperature estimates in the absence of long-term meteorological stations (Jacobs and others 2014; Finnis and Bell 2015; Way and Bonnaventure 2015; Way and Viau 2015; Way and others 2017). Daily mean surface air temperatures (SATs) from ERA5 atmospheric reanalysis (Hersbach and others 2020) were extracted using Google Earth Engine (Gorelick and others 2017) for grid cells (0.25°) covering most of TMNP and adjacent areas (January 1979–August 2019). ERA5-derived SATs were then locally scaled with linear regression using short records of SAT from ecosystem monitoring stations operated by Parks Canada and Centre d’études nordiques (Sarrazin and Allard 2020) in northern Labrador and Nunavik (N = 15; mean coverage of 5.7 years per station since 2010). Daily SATs for Nakvak Brook (427 m asl) were then calculated from the locally scaled ERA5 estimates using thin-plate spline interpolation within the ‘fields’ package (Nychka and others 2017) in R (R Core Team 2020), employing positional coordinates and elevation as predictors (as in Way and others 2017).

Time series of the Normalized Difference Vegetation Index (NDVI) derived from remote-sensing imagery are commonly used as a proxy of decadal-scale changes in vegetation productivity. In high-latitude environments, NDVI is related to above-ground phytomass (Raynolds and others 2012) with shrubs generally having higher NDVI values due to their greater canopy leaf area relative to other tundra plants (Boelman and others 2011). Landsat imagery (5, 7, and 8; 30 m resolution) was used to develop an NDVI time series for a 250 km2 area around the study site at Nakvak Brook (Figure 1B). A quality control procedure developed by Pironkova and others (2018) automated the acquisition of cloud-free, sensor-calibrated images in Google Earth Engine from which annual median July–August NDVI values were calculated. Due to persistent cloud cover in some years, site-level median NDVI values are available for only 23 of 34 years. Estimates of monotonic rates of change over 1985–2019 were calculated using Theil–Sen’s slope estimator in the ‘EcoGenetics’ package (Roser and others 2017) in R. The distributions of NDVI values between averages of the first (1985–1994) and most recent (2010–2019) 10 years of observation were compared using a Kolmogorov–Smirnov test. Together, the local temperature and NDVI time series provide a decadal-scale record of environmental change that contextualizes recent field-based observations.

Plot-Scale Changes in Plant Growth and Composition

A series of vegetation surveys were undertaken to connect regional-scale changes with observations at the local level. Surveys of shrub composition were performed along three transects at Nakvak Brook in the summer of 2017 (see site and transect layout in Figure 1C). Shrub ramets were collected from a representative sample of individuals in each transect (n = 4–12) to develop establishment and growth records. Care was taken to ensure that ramets were collected from unique plants because many arctic shrubs reproduce vegetatively through the production of rhizomes. The largest, and presumed oldest, shrub ramet was cut at the root collar, and the basal cross sections were used to estimate the minimum year of establishment and to develop an annually resolved ring-width chronology. Samples from Transect 1 were also serially sectioned at 10-cm intervals along the remainder of the stem to calculate rates of stem elongation.

Ramet cross sections were prepared, dated, and analysed using standard dendrochronological methods adapted for shrub species (Myers-Smith and others 2015; see Supplemental Material 1). Recent increases in shrub size were determined from ramet subsections (Figure 1C; Transect 1) by examining changes in the rates of stem elongation over time. Stem elongation data were used to fit linear and negative exponential models with the position of the sample along the ramet (cm from base) as the dependent variable and the year in which it reached that length as the independent variable. A value of 0.1 cm was added to stem lengths in the exponential model to permit the inclusion of the initial 0 cm class in the model. The model with the highest explained variance (greatest R2) was considered to have the best fit, with a linear model implying a constant rate of stem elongation over time and a negative exponential implying an increasing rate. Finally, an annual radial growth chronology was established for the study area using the growth measurements of the basal samples of B. glandulosa, the most common shrub species in the transects. The chronology was developed from an initial 18 individuals by detrending cross-dated ring-width series using the series mean, as no age-related growth trends were apparent (Myers-Smith and others 2015). Detrended series were combined using Tukey’s bi-weight robust mean, and autoregressive pre-whitening was applied to account for non-independence between subsequent years of shrub growth (Cook and others 1990). Spearman rank correlations were performed (due to limited years of observation) between the residual B. glandulosa growth chronology (1998–2017; 20 years) and three temperature records (mean July, August, and July–August temperatures); between July–August NDVI and the growth chronology (1998–2017; 14 years due to missing years of NDVI data); and between NDVI and the temperature records (1985–2017; 23 years due to missing years of NDVI data).

Point-frame vegetation surveys were conducted during the 2010 and 2015 growing seasons to identify short-term temporal changes in plant composition and abundance at Nakvak Brook under ambient climate conditions. Methods developed for the International Tundra Experiment (ITEX; https://ibis.geog.ubc.ca/itex/) were used to collect information about plant cover, size, and abundance in permanent sample plots. To account for the moderating effect of soil moisture on plant responses to climate warming (Elmendorf and others 2012b; Ackerman and others 2017; Bjorkman and others 2018), a representative sample of 10 ‘wet’ and 10 ‘dry’ plots were established in 2007/2008. Wet plots were qualitatively distinguished as having saturated soil and/or some standing water, whereas dry plots had no standing water and had dry soils. Instantaneous soil moisture recordings were made during field visits in July or August of 2009, 2010, 2011, and 2013 using a type wet-2 WET Sensor and HH2 Moisture Meter (Delta-T Devices Ltd) (Supplementary Figure 1). The plots measured 1 m × 1 m and the point-frame had a grid spacing of 10 cm, resulting in 100 points for vegetation data collection per plot. During the surveys, each encounter between a plant and a point-frame pin, and the height at which it occurred, was recorded. Plants were identified to the species level when possible; however, some could only be identified to genus. Plants were grouped as deciduous shrub, evergreen shrub, forb, graminoid, lichen, or moss for statistical analyses. Bare ground was categorized as rock, exposed soil, or cryptogamic crust, and standing dead plants or litter were grouped as dead organic material.

Linear mixed-effect (LME) and generalized linear mixed models (GLMM; ‘lme4’ package; Bates and others 2015) were used to evaluate changes in bare ground, plot canopy height, maximum plant height, living plant abundance, and abundance of dead plant material across years and between plot moisture status (see Supplemental Material 2). LME models were used to fit maximum plant (one model per lifeform) and mean canopy height data, and GLMMs were used for bare ground cover (binomial distribution), life-form abundance (one model per lifeform; Poisson distribution with log-link), and dead plant material data (Poisson distribution with log-link). In all models, soil moisture class (wet versus dry), survey year (2010 versus 2015), and their interaction were the fixed effects, and a categorical variable for plot ID was the random term.

Ground surface temperatures (GST) were monitored at the corners of vegetation plots from 2010 to 2019 using HOBO(R) Onset UA-001-08 pendant temperature loggers to infer changes in local surface conditions over the observational period and between wet and dry plots. The temperature loggers were buried between 2 and 5 cm below the soil surface and recorded temperatures at 1-h or 2-h intervals, depending on the year. Ground surface data were interpolated to an hourly interval (maximum gap of 3 h) using a cubic spline to correct for time of observation bias and were then aggregated to daily mean values. Ground surface temperature data were analysed using the temperature at the top of the permafrost model framework (Smith and Riseborough 1996; Riseborough 2004) following Way and Lewkowicz (2018). Interannual differences between warm season air (local ERA5 estimates) and soil (observed) temperatures were summarized using the thawing season offset (TSO; for example, Way and Lewkowicz 2018), whereas differences between cold season air (estimated) and soil (observed) temperatures used the nival offset (NVO; for example, Riseborough 2004; Way and Lewkowicz 2018 see Supplementary Material 3). Using TSO and NVO is preferred to more commonly used metrics (for example, n-factors) because it enables direct comparisons of net air-surface temperature differences during warm and cold seasons, respectively. Statistical differences in TSO and NVO between wet and dry plots were determined for 2012–2018 (years with sufficient data coverage) using linear mixed-effects models; soil moisture class, year, and their interaction were fixed effects, and plot ID was the random term. A Spearman correlation was performed to determine if relationships existed between average annual NVO, TSO, and the average height of the plot canopy.

Permafrost Investigation

Permafrost characteristics at Nakvak Brook were investigated in late July 2017 using DC electrical resistivity tomography (ERT) and standard field methods including frost probing, instantaneous temperature profiling, and ground temperature observations (Lewkowicz and others 2011). ERT can be used to infer the presence or absence of frozen soil or rock by modelling the pattern of measured resistance to electrical current between pairs of electrodes inserted in the ground (Hauck and others 2003; Hauck 2013). Frozen ground is generally more resistive to current than unfrozen materials due to lower unfrozen moisture contents (Hauck 2013), but dry unfrozen soils and bedrock can exhibit resistivities that overlap with the range observed for frozen soils. Given the relatively coarse soils and near-surface bedrock present at Nakvak Brook, validation in relation to frozen ground conditions proved challenging.

ERT investigations were conducted along three transects that sampled a variety of local conditions (Figure 1C) and were co-located with shrub surveys. Resistivities were measured with an ABEM Terrameter LS profiling system using a Wenner configuration across 40 m (N = 2) and 160 m (N = 1) transects (maximum depth penetration of 7 m and 26 m, respectively). RES2DINV was used to invert measured resistivities with the robust inversion method (Loke and Barker 1996, 2003) and was subsequently exported to xyz data after errors dropped below 5% and/or error values fell by < 1% between model fit iterations (for example, Way and others 2018). Profiles were topographically corrected using elevations collected in the field with a handheld GPS (Garmin Oregon 450t) and a Brunton compass. ERT tomograms are presented as model blocks with larger blocks at depth and near the ends of profiles. Tomograms were generated using a customized R v3.4 script, and areas below the depth of investigation from RES2DINV (per Loke and others 2003) were excluded. For validation, we used instantaneous ground temperature profiles collected with Onset Hobo UX120-006M Analog Data Loggers (accuracy ± 0.15 °C) connected to two to four vertically arranged thermistors separated by 20 cm and mounted on a wooden dowel (for example, Way and Lewkowicz 2015; Holloway and Lewkowicz 2020). The observed thermal gradient in the lower 20 cm was used to extrapolate temperatures to 0 °C in order to approximate the thaw depth. This likely underestimates the active layer thickness because measurements were made before the end of the thaw season.

To further elucidate permafrost conditions at Nakvak Brook, 1000 simulations of permafrost probability were run using the equilibrium temperature at the top of the permafrost model (Smith and Riseborough 1996) for 14 vegetation plots that had longer records of ground surface temperatures (see Supplementary Material 2 for additional details regarding permafrost simulations).

Results

Decadal-Scale Changes in Air Temperature and Vegetation

Locally scaled ERA5 data showed warming trends at Nakvak Brook in all seasons (Figure 2). Cooling was observed from 1979 to 1990, followed by rapid warming until 2011, and then elevated but variable temperatures from 2011 to 2019. Temperature increases in autumn and winter were greater than those in spring and summer.

Figure 2
figure 2

Changes in seasonal average temperatures for Nakvak Brook (1979–2019) derived from locally scaled ERA5 atmospheric reanalysis data. Ordinary least squares regression lines are shown with shaded areas indicating 95% confidence intervals (note different temperature ranges on y-axes).

The median July–August NDVI record derived from Landsat showed a general increase over time. Rates of change in NDVI (Theil–Sen slope estimates; change in NDVI value per year) varied spatially, with the greatest increase occurring along the valley-side slopes of drainage channels now occupied by tall, upright shrubs (south-west portion of Figure 3). Average NDVI values increased from 0.218 to 0.328 between the first (1985–1994) and most recent (2010–2019) ten-year periods of observation, and the Kolmogorov–Smirnov test indicated a significant difference in cumulative distribution functions (D = 0.260; p < 0.001) as values became more positive over time.

Figure 3
figure 3

Observed changes in NDVI from 1985 to 2019 shown for a 250 km2 area around the study site at Nakvak Brook (red box, centre). Rates of change are significant Theil–Sen slope estimates (a < 0.05) measured in NDVI units per year. Bottom right: Frequency histogram of the distribution of NDVI values during the first (1985–1994; light green) and most recent (2010–2019; dark green) 10 years of observation, for years in which adequate imagery was available. Dashed coloured lines represent the mean NDVI value of each time period.

Plot-Scale Changes in Vegetation Growth and Composition

Shrubs growing in the survey transects were primarily dwarf birch (B. glandulosa, 18 of 26 individuals), with the remainder being species of willow (Salix pedicellaris, 7 of 26; Salix sp., 1 of 26). Ramet establishment varied between 1997 and 2014 and peaked in 2006 (Supplementary Figure 3a). The oldest shrubs were B. glandulosa and all Salix established after 2005. Shrub establishment tended to be greater in years with higher average July temperatures, although this observation is based on a small number of samples (Supplementary Figure 3b). Serial sectioning of shrub ramets from Transect 1 revealed that stem elongation increased exponentially over time at the transect-level (R2 of negative exponential model = 0.64; p < 0.001 relative to R2 = 0.50; p < 0.001 for the linear model; Figure 4). Modelled shrub stem lengths increased rapidly after ca. 2005, with annual growth rates increasing from less than 1.00 cm/y in 2005 to 13.25 cm/y in 2014.

Figure 4
figure 4

Rate of stem elongation (cm/yr) for 10 individuals growing in Transect 1 with colours indicating shrub species. Elongation rates were best described by an exponential relationship, indicating that the transect-level growth rate of dominant stems has accelerated over time. A slight jitter in points along the x-axis was used to show overlapping data.

The final radial growth chronology contained measurements from 13 dwarf birch ramets and spanned 21 years (1997–2017; Supplementary Figure 3b). Shrub radial growth generally increased and became more variable over time. Although performed on a limited number of data points (14–23 years), correlation analyses between shrub growth, NDVI, and summer temperatures show significant, positive relationships (Figure 5). The strongest correlations were with mean August air temperatures, rather than mean July or July–August values (p > 0.05). August is typically the warmest month of the growing season at Nakvak Brook, as noted above.

Figure 5
figure 5

Relationships between: A shrub ring widths and average August temperatures; B shrub ring widths and median July–August NDVI; and C July–August NDVI and August temperatures for the study area (*p ≤ 0.05; **p ≤ 0.01).

Interestingly, the point-frame vegetation surveys showed plant size and abundance to differ most between moisture categories (wet versus dry plots) as opposed to between survey years (2010 and 2015). Plant canopies were taller in wet than in dry plots (twet = 4.61; p < 0.001; Figure 6; see Supplementary Tables 1–3 for full model parameter estimates), as were the maximum heights of deciduous shrubs (twet = 10.70; p < 0.001), forbs (twet = 6.50, p < 0.001), and graminoids (twet = 8.60; p < 0.001). Forbs (zwet = 3.13; p < 0.001), graminoids (zwet = 2.00; p < 0.001), and moss (zwet = 1.47; p < 0.001) were more abundant in wet plots, whereas evergreen shrubs and lichen were favoured in dry plots where they were both taller (evergreen shrubs: twet = −3.20; p = 0.03) and more abundant (evergreen shrubs: zwet = −4.18; p < 0.001; lichen: zwet = −4.62; p < 0.001).

Figure 6
figure 6

Effects of moisture class, year, and their interaction on: A the maximum height of each life-form; B life-form abundance; C average height of the plot canopy; D abundance of dead plant material; and E bare ground cover. Significant effects have 95% confidence intervals (bars) that do not overlap zero. Dry plots and the 2010 survey are the baseline categories against which the effects are measured.

Temporal changes in plot surveys mainly relate to differences in life-form abundance, although the maximum height of graminoids uniquely decreased from 2010 to 2015 regardless of soil moisture status (t2015 = −3.10; p = 0.010). Forbs increased in abundance in 2015 (z2015 = 1.27; p = 0.003), and more so in dry than wet plots (zwet*year = −1.09; p = 0.015). This change in abundance was driven by an increase in the amount of some species (for example, Polygonum viviparum) as well as the appearance of new species that were not present in the plots in 2010 (for example, Equisetum arvense, Euphrasia arctica; Oxytropus campestris). In contrast, the interaction between plot moisture category and year was significantly positive in the model of graminoid abundance (zwet*year = 0.43; p = 0.013), indicating that the increase was greater in wet plots. Lichen abundance increased very slightly from 2010 to 2015, independent of moisture class (z2015 = 0.27; p = 0.001), and moss declined in wet plots in 2015 (zwet*year = −0.59; p < 0.001). The likelihood of encountering standing dead plant or plant litter decreased in 2015 relative to 2010 (zyear = −0.88; p < 0.001), and declined more in wet than dry plots (zwet*year = 0.44; p = 0.001). Finally, the GLMM for bare ground cover showed there to be no differences between years (zwet = 0.18; p = 0.282) or between plot status (zwet = −1.22; p = 0.170). There was no evidence of caribou herbivory within the plots or in the immediate area.

Ground surface temperatures recorded at vegetation plots showed that air-surface temperature differences were larger in the winter (that is, higher NVO) at wet sites than at dry sites (Figure 7; Supplementary Table 4). Snow cover warmed the ground surface by about 2–6 °C relative to air temperatures, and the average difference in NVO between wet and dry plots was 1.15 °C (LME; pwet = 0.011). Wet plots showed less variability in NVO than dry plots, indicating greater consistency in the air-surface temperature decoupling in the cold season at wet sites. The correlation analyses showed mean NVO to be positively related to plot canopy height (Rho = 0.68; p = 0.010). Unlike NVO, TSO did not differ significantly between moisture classes (twet = 0.04; p = 0.805) and its impact on GSTs was smaller, generally between −0.5 °C and + 0.5 °C. This indicates that shading by vegetation results in minimal air-surface temperature decoupling in both wet and dry plots, with thawing season air temperatures warmer than the ground surface in some years and colder in others, and with no consistent trend. The relationship between TSO and canopy height was negative but not significant (Rho = −0.37; p = 0.200). Overall, the influence of NVO on GST, mainly due to snow, was often an order of magnitude larger than the influence of TSO, mainly due to shading by plants and other thawing season impacts (for example, soil moisture). Thus, interannual variability in the yearly offset between air and GST is primarily driven by differences in snow accumulation and redistribution around the study site.

Figure 7
figure 7

Annual differences in A nival (NVO) and B thawing (TSO) season offset between wet and dry vegetation plots (2012–2018). NVO and TSO values were derived from ground surface temperature loggers and locally scaled air temperatures for Nakvak Brook. Differences between wet and dry plots were determined using Tukey post hoc contrasts (*p ≤ 0.05; **p ≤ 0.01). Note differences in scale of y-axis.

Site-Specific Permafrost Investigations with ERT

Transect 1

Transect 1 traversed a very gently sloping section of the valley bottom covered by low herbaceous vegetation with prostrate willows in places. Small (10–20 cm) rises and falls in the surface corresponded to locally drier and wetter conditions in the fine-grained surficial deposits. Frost-shattered rocks on the surface near the start of the transect suggest that bedrock contact may be quite close to the surface. The tomogram shows a 0.5–1-m-thick layer of mostly low resistivities (minimum: 1075 Ωm; median: 2100 Ωm) overlying a body of higher resistivities (minimum: 1840 Ωm; median: 13,700 Ωm) that extends to the base of the profile at 7 m depth (Figure 8A). A clear frost table was encountered at depths ranging from 28 to 50 cm at several locations along the transect during probing in late July. Instantaneous ground temperatures collected at 12 m and 35 m along the profile also indicate the presence of frozen ground at depths of 47 cm and 40 cm. A permafrost mound about 70 cm high with a thaw depth estimated at 43 cm was present within 35 m of the profile. Taken together, the observations at Transect 1 indicate a thin, thawed layer overlying frozen ground, some of which could be bedrock. Given the timing of the survey and the shallow depths of thaw, these are interpreted as representing an active layer of variable thickness overlying permafrost at least 7 m thick.

Figure 8
figure 8

XYZ plots of modelled resistivities acquired with ERT at A Transect 1, B Transect 2 and C Transect 3 in late July 2017. Interpreted versions of the ERT profiles are presented in Supplementary Figure 4. Note elevation coordinates in B are not available.

Transect 2

Transect 2 ran across the valley bottom, and in the first 20 m and last 5 m, it intersected a series of drainage channels. The channels were mostly < 1 m wide, incised 30–40 cm relative to the adjacent terrain and lined with angular rocks 15–30 cm in diameter. Water was present in the channels, and soil in the areas between them was saturated at a depth of 30 cm. Vegetation along the transect was classified as wet sedge. The tomogram shows a two-layered system, with a 2–3-m-thick low resistivity layer (minimum: 112 Ωm; median: 745 Ωm) with some higher resistivity patches linked to problems with contact resistance, bounded by a sharp gradient with a higher resistivity body (minimum: 2010 Ωm; median: 18,520 Ωm) that extends to the base of the profile at 7 m depth (Figure 8B). Probing for the frost table did not reveal a defined thaw front along this profile, as clasts were encountered at shallow depths in many places. The maximum probe depth reached was 88 cm at 34 m along the profile, but refusal was due to a clast rather than frozen ground as the temperature at this depth was 7 °C. Two instantaneous ground temperature measurements gave estimated depths for 0 °C as 251 cm and 384 cm. Considering the substrate and drainage conditions at Transect 2, we infer that the low resistivity values in the near-surface indicate soil thawed to at least a depth of 2–3 m, and saturated conditions. High resistivities at greater depths could either be deep permafrost or bedrock. Even if permafrost is present, the depth of the 0 °C isotherm attained by mid-summer, as indicated by the tomogram and the ground temperature measurements, shows it is likely that the thawed layer would not completely refreeze in winter, meaning a talik would be present.

Transect 3

The third transect ran up the foot of a WNW-facing slope that rises to an elevation of about 800 m asl over a distance of 1.9 km. Only the basal 160 m of this slope was surveyed. Surficial materials along the transect comprised angular rocks embedded in a matrix of coarse sediments, interpreted as till and colluvium. The first 110 m of the profile was cut by a network of shallow channels lined with rocks as at Transect 2, and with water levels 25–30 cm below the adjacent terrain. Between 110 and 160 m, the slope became drier and slightly steeper. Vegetation was comprised of grasses and mosses in the lowest part of the profile with scattered willows. Cover decreased in the last section of the profile which was drier and had a higher concentration of surface rocks.

The tomogram for Transect 3 is the most complex of the three. In general, resistivities increase with depth, but there is considerable lateral variation in the near-surface layers (Figure 8C). Values in the first half of the profile are as low as 1300 Ωm in the upper 3 m close to or in the channels, rising to > 13,000 Ωm at the base of the profile at 25 m depth. The near-surface to depth resistivity gradient is much reduced in the upper half of the profile where the slope gradient is about 5° and surficial cover appears to become thinner. Near-surface resistivities range from 2150 to 9550 Ωm, while deeper layers are in the 11,000 Ωm–15,000 Ωm range. Frost tables could be detected by probing at only 8 of 40 locations along the profile due to the high concentration of clasts, and all these were in the first half of the transect. Thaw depths ranged from 57 to 95 cm, similar to two values (75–90 cm) estimated from instantaneous ground temperatures measured in the same parts of the profile. Probing depths attained in the final 40 m of the transect were only 10–40 cm which demonstrates high clast contents.

Combining all the data collected for this transect, it is concluded that permafrost is likely present throughout, but surficial sediments become thicker downslope which allows the wet surface thawed layer to be distinguished in the resistivity profile. The gradational increase in resistivities to a depth of 10 m in this part of the profile suggests changes in frozen ground conditions (possibly due to temperature or ice content change). Below 10 m depth, the high resistivities could be due to more frozen soil layers, or bedrock. In the upper part of the profile, drier conditions and from the geomorphological context, a possible bedrock contact approaching the surface, resulted in higher resistivities and less contrast with depth.

Simulation of Permafrost Probability

Simulated permafrost probabilities based on ground surface temperature data at the vegetation plots ranged from 9 to 100% with a mean probability of permafrost of 55% for the 14 logger sites (Figure 9). Nine of 14 sites had permafrost probabilities exceeding 25%, and dry sites typically had a much higher permafrost probability and less variability (mean: 79%; sd: 30%) compared to wet sites (mean: 32%; sd: 28%). These data indicate the distribution of permafrost at Nakvak Brook is widespread but discontinuous, in agreement with our interpretations of the ERT data.

Figure 9
figure 9

Simulated permafrost probability for 14 of 20 vegetation plots estimated from 1000 simulations using the temperature at the top of the permafrost model. Models are based on mean ground surface temperatures over the period 2011–2018.

Discussion

Using a mixed-methods approach, we found evidence of rapid, ongoing changes in climate and vegetation conditions in the low Arctic Torngat Mountains of northern Labrador. We have attempted to move towards an integrated understanding of interactions between above-ground and subsurface conditions, which are complex and difficult to measure directly. Further environmental changes are expected in the coming decades, and our results highlight how and where future changes may manifest in the region.

Relationships Between NDVI, Climate Change, and Shrub Growth

Observed increases in NDVI at Nakvak Brook align with other reports of rapid vegetation change in the eastern Canadian Arctic and Subarctic (Fraser and others 2011; Tremblay and others 2012; Quirouette and Zorn 2015) and demonstrate that greening has continued since the publication of these earlier studies. Tundra greening is generally related to increased plant productivity or changes in species composition (Myers-Smith and others 2020); our dendrochronological analyses and point-frame surveys provide insight into the factors underlying the increase in NDVI at Nakvak Brook.

Recent shrubification is indicated by the records of ramet establishment and shrub growth. The ramet establishment dates show that stem initiation began in 1997 and peaked in 2006. Arctic shrubs have relatively long potential lifespans (decades to over a century; Myers-Smith and others 2015 and references therein); thus, the young ages of the ramets could indicate that limiting conditions were relaxed relatively recently. We did not quantitatively analyse the relationship between climate and establishment frequency; however, stem initiation tended to be higher during years with warmer July temperatures suggesting that a pulse in shrub establishment may have been triggered by warming in the mid-to-late 1990s. Shrub recruitment elsewhere has shown variable responses to climate. Myers-Smith and Hik (2017) found recruitment to be positively correlated to winter temperatures (willows; Yukon, Canada), whereas Frost and Epstein (2014) found tall shrub expansion to be more sensitive to precipitation (alder; northern Siberia). Büntgen and others (2015) found historical shrub recruitment to initially respond to increased temperatures but to decline in recent decades despite renewed warming (various species; East Greenland). Variability in recruitment responses to climate change may be caused by regional and species-level differences in climate limitations and the moderating influence of local factors (for example, herbivory; soil moisture; disturbance; Frost and others 2013; Lantz and others 2013). A comprehensive study of the relationships between climate and shrub recruitment would help confirm whether summer temperatures are a dominant control on shrub expansion in the study region.

Increases in NDVI at Nakvak Brook also appear to be driven by recent and rapid shrub growth. The lengths of shrub ramets increased exponentially at the site level, indicating an easing of limitations to primary growth over time. Stem lengths of Arctic shrubs have been shown to respond favourably to ambient warming and to experimental warming and fertilization treatments (Bret-Harte and others 2002). The analysis of our dwarf birch growth chronology also demonstrates a quantitative link between shrub growth, summer temperatures, and NDVI. Ropars and others (2015a) similarly found the radial growth of dwarf birch to be positively correlated to summer temperatures and NDVI in Nunavik, and suggested the greening trend detected in that region to be driven by dwarf birch growth.

The focus of our research was on climatic impacts; however, it is possible that reduced herbivore pressure from caribou in the Torngat Mountains contributed to the process of shrubification as the large decline in numbers between the 1980s and 2014 coincided with a period of warming. Research has shown that herbivory by Arctic ungulates can reduce but generally not reverse the process of tundra shrubification (Christie and others 2015). For instance, in Nunavik, Morrissette-Boileau and others (2018) found that herbivory did not counteract the influence of climate warming on shrub expansion, even during periods of very high caribou abundance. Experimental field studies and continued caribou monitoring are required to determine if vegetation in the Torngat Mountains is similarly resilient to the impacts of caribou browsing.

Plot-Scale Vegetation Change

Changes in vegetation detected in point-frame surveys were more subtle than the temporal trends observed in the shrub-ring records, and depended heavily on site context. Plot-level differences in vegetation size and abundance were mainly driven by the moisture status of the plots, reflecting the importance of microhabitat conditions in determining broader patterns of vegetation composition (Oberbauer and Dawson 1992). Previous research has found soil moisture to have a moderating effect on the sensitivity of tundra plant responses to temperature warming across multiple spatial scales (Elmendorf and others 2012a; Ackerman and others 2017; Bjorkman and others 2018), and we also found some indication of this tendency at Nakvak Brook.

Interactions between soil moisture and time resulted in a decline in moss abundance in 2015 in wet plots; mosses are sensitive to soil moisture and their tolerance to desiccation can be reduced by prolonged exposure to warm and dry conditions (Dorrepaal and others 2003). Our field observations of soil moisture suggest a drying trend among wet plots at Nakvak Brook, which combined with elevated temperatures, likely contributed to the decline in mosses. A surprising finding of the point-frame surveys was a slight but statistically significant increase in lichen abundance between years that was not related to plot moisture (although lichens were more abundant in dry plots overall). Vegetation surveys have frequently reported reduced lichen cover in response to warming temperatures (Cornelissen and others 2001; Walker and others 2006; Paradis and others 2016). Similar to our findings, however, Siegwart Collier (2020) observed lichen cover to increase from 2009 to 2016 in her vegetation plots exposed to ambient conditions in Nain and Torr Bay, Nunatsiavut. Lichen decline is typically attributed to increased shading by deciduous shrubs and plant litter, which we did not observe at Nakvak Brook (that is, deciduous shrubs were not significantly taller or abundant in 2015 than in 2010, and litter/dead plant material declined).

Altogether, the effect of site moisture class on vegetation cover and abundance was stronger than the effect of sampling year, which may be due to the lack of a warming trend over the 5 years between observations. The strength of warming has been shown to influence the likelihood of plants undergoing an increase in abundance or cover (Elmendorf and others 2012b). Locally scaled climate data showed summer temperatures to be higher in 2010, the first year of data collection, than in any year since 1979. Although remaining elevated above the long-term average, temperatures underwent a short-term decline from 2010 to 2015.

Connections Between Vegetation Change and Subsurface Characteristics

Ground surface temperature data provide a link between atmospheric temperatures, plant growth, and the subsurface conditions of the study site. Plot-level comparisons of NVO and TSO between moisture classes show that wet plots, in which deciduous shrubs were taller, experienced higher cold season soil temperatures than dry plots. The ability of shrubs to trap snow and insulate the ground has been widely reported (Pomeroy and others 1997; Sturm and others 2005; Myers-Smith and Hik 2013; Roy-Léveillée and others 2014; Domine and others 2016), and we also found the insulating effect to be moderated by soil moisture, possibly through its influence on shrub size or soil freezing characteristics. Annual NVO values at Nakvak Brook varied by an order of magnitude (0.8 °C–8.4 °C) and were comparable to other reported winter warming effects (Sturm and others 2001a, 2005; Myers-Smith and Hik 2013; Domine and others 2016; Paradis and others 2016). The more modest insulating effect in some plots at our study site could be due to the low overall size of the shrubs, since the insulating effect of snow has been shown to increase with plant height (Roy-Léveillée and others 2014). Alternatively, the greater maximum cooling effect in wet soils due to the release of latent heat during freezing may have contributed to differences in ground surface temperatures between plot types (Romanovsky and Osterkamp 2000). In either case, plot canopy height was positively correlated to average NVO values, and we expect that the insulating effect of shrub cover may increase in the future if shrubs continue to increase in size.

Analysis of in situ temperature and ERT data combined with field observations of permafrost (palsas) and frost-associated phenomena (sorted and polygonal features) suggest that permafrost is widespread at Nakvak Brook (Supplementary Figure 5), but generally at depths below the rooting zone of shrubs. However, the presence of near-surface permafrost (for example, upper 3–4 m) is challenging to interpret within the study area. At Transect 1, inferred permafrost was ubiquitous along the ERT profile, whereas at Transect 2 it was inferred to be absent. At Transect 3, frozen ground in the near-surface layers was likely present along sections of the transect, while other sections were difficult to interpret or impacted by near-surface bedrock. Stream channels were areas of localized deeper thaw and permafrost was largely absent in the near-surface. Generally, coarser surficial materials and very wet surface conditions were associated with an inferred deeper active layer and warmer ground temperatures in the near-surface layers while finer surface materials were associated with a thinner active layer and moderate surface wetness.

Field observations of surface wetness changes during annual field visits to Transect 2 in particular suggest that subsurface hydrological changes may have occurred prior to enhanced shrub invasion. The contrast between frozen ground with a shallow active layer (Transect 1) and absent near-surface permafrost (Transect 2) makes it plausible that permafrost thaw could have modified subsurface hydrological conduits at the latter site over the period of observation. We hypothesize that the potential for ecological impacts of permafrost change at Nakvak Brook depends heavily on geomorphological context with moderately wet sites having finer materials being the most likely to have permafrost currently limiting vegetation growth and representing potential hot-spots for future change. In general, the permafrost simulations combined with the ERT data suggest that mid-elevations in the southern portions of the Torngat Mountains may be better classified as falling within the extensive discontinuous permafrost zone per Way and Lewkowicz (2016) as opposed to the continuous permafrost zone (for example, Permafrost Map of Canada).

Implications of Change

Over the past four decades, significant changes have occurred at Nakvak Brook with implications both for people as well as the broader environment. The changes we have documented may affect how Nunatsiavummiut and Nunavimmiut use the region for travel, hunting, and gathering plants. Shrubification in Torngat Mountains National Park will make it more difficult for Inuit to access travel routes by snowmobile, may alter the distribution of berry plants, and will change food availability for caribou and bear. In valley bottoms, where we observed significant greening, tall shrubs also provide cover for black bears and polar bears, increasing safety concerns. The study area is also located within one of Canada’s National Parks and is intended to be representative of the southern limit of the Canadian Arctic Cordillera. Many protected areas in North America are expected to experience a shift in biome representation by the end of the century (Lemieux and Scott 2005; Holsinger and others 2019), and northern protected areas may become increasingly important refuges for species migrating northward (Berteaux and others 2018). Our results indicate that environmental change is already detectable at an extreme, mid-elevation research basin in Torngat Mountains National Park, an example of the disproportionate influence of climate change on protected areas (for example, Gonzalez and others 2018).

Conclusions

Applying an interdisciplinary lens and mixed-methods approach, this study provides a multivariable assessment of environmental change at the Nakvak Brook research basin located in the low shrub tundra ecotone in Torngat Mountains National Park, northern Labrador. Our findings confirm that vegetation greening is ongoing in the area and provide direct evidence that greening is driven by the influence of warming temperatures on shrub establishment and growth. Our data establish a link between vegetation cover and ground surface temperatures during the snow season with the expectation that future growth responses to warming may enhance this effect. At the plot scale, we found plant size and cover differences to be primarily driven by soil moisture regime, a characteristic that is tied both to atmospheric warming as well as subsurface conditions and that will be an important mediator of future responses to climate change. Near-surface permafrost was variable throughout the study area, and its importance for above-ground vegetation depends heavily on local geomorphological context. Continued research into connections above and below the soil surface and regional hydrological connectivity is warranted to evaluate the impacts of climate change on vegetation and permafrost conditions over the coming decades.