Study site

The study site is located in a 2 km² watershed in the Tasiapik Valley near the Inuit community of Umiujaq along the east coast of Hudson Bay in Nunavik (Québec), Canada (Figure 2), which lies within the discontinuous permafrost zone (Allard and Seguin 1987). Following the retreat of the Laurentide Ice Sheet and marine transgression around 8000 cal years BP (Hilaire-Marcel 1976; Lavoie et al. 2012), different types of Quaternary sediments were deposited in the valley floor such as coarse-grained glacial and fluvio-glacial sediments, fine-grained offshore sediments, and intertidal and littoral sediments (Fortier et al. 2020). Due to isostatic rebound after deglaciation, the Quaternary deposits in the valley gradually emerged and came in contact with the subarctic cold climate, which induced permafrost aggradation and formation of permafrost mounds (Figures 1 and 2c) (Lafortune et al. 2006; Lavoie et al. 2012; Fortier et al. 2017, 2020). Evidence of permafrost degradation including thaw settlement, formation of thermokarst ponds, pellicular slope movement, and shrubification is already seen in the Tasiapik Valley due to climate warming.

Climate

The time series of mean annual air temperature (MAAT) at Umiujaq from 1926 to 2019 is provided in Figure 3. Because air temperature in the Tasiapik Valley has only been monitored since 2000 at a height of 3 m above ground at a meteorological station known as VDT-SYBU, the time series was completed based on the climate variability recorded at meteorological stations operated by Environment Canada in the neighboring Inuit communities of Kuujjuarapik and Inukjuak, located along the east coast of Hudson Bay, south and north of Umiujaq, respectively. The MAAT at Umiujaq (MAAT_UMIUJAQ) from 1926 to 2000 was assessed by averaging the MAATs at Kuujjuarapik and Inukjuak (MAAT_KUUJJUARAPIK and MAAT_INUKJUAK). This climate variability at Umiujaq (Figure 3), which controls permafrost dynamics, is considered in the numerical simulations presented and discussed in the next sections.

Figure 2.

(a) Location of the Inuit community of Umiujaq, Nunavik, Canada. (b) IKONOS satellite image showing the Tasiapik Valley and study site location. (c) Permafrost mounds 1 and 2 at the study site.

Outline

The numerical model was designed to allow ground surface conditions to evolve in space and time as the climate warms. As detailed below, ground surface temperatures corresponding to these conditions were derived from simple linear regressions between measured ground surface temperatures and air temperatures over a full year for different characteristic ground surface conditions. Assuming these relationships remain constant in the future as the climate warms, hypothetical scenarios of future increases in air temperature are then used to predict corresponding future changes in ground temperatures due not only to climate warming but also to changes in ground surface conditions which drive the numerical model and permafrost dynamics.

Characteristic ground surface conditions and associated ground surface temperature variability were assessed from infrared thermal imaging of permafrost mounds 1 and 2 (Figure 2c) using a FLIR T400 infrared camera. Examples of infrared photographs of the study site are provided as supplementary material ( Figures SM1 (teco_a_1949819_sm2254.zip)).

Surface temperature monitoring

Onset HOBO™ Pro v2 temperature probes were buried 5 cm below the ground surface to monitor surface temperatures of each characteristic ground surface condition observed in the field and delineated with infrared thermal imaging ( Figure SM1 (teco_a_1949819_sm2254.zip)). The accuracy of these temperature probes is ±0.2°C from 0°C to 50°C while the resolution is 0.02°C at 25°C. In total, 25 probes were installed at the study site ( Figures SM2 to SM4 (teco_a_1949819_sm2254.zip)) which recorded hourly temperatures from July 2018 to July 2019.

Four examples of time series of hourly air and ground surface temperatures, and graphs of mean daily ground surface temperature as a function of mean daily air temperature, are provided in Figures 4 and 5 for different ground surface conditions: (1) mosses and lichens (temperature probe TR-6; Figure 4a and 4c), (2) shrub vegetation (temperature probe TR-1; Figure 4b and 4d), (3) initial stage of a developing thermokarst pond (temperature probe MT-1; Figure 5a and 5c), and (4) final stage of a thermokarst pond (temperature probe MT-3; Figure 5b and 5d). Temperature probe TR-6 was covered by 0.1 m high mosses and lichens ( Figure SM4a (teco_a_1949819_sm2254.zip)) whereas temperature probe TR-1 was covered by dense and tall shrubs up to 2 m high ( Figure SM4a (teco_a_1949819_sm2254.zip)). Temperature probes MT-1 and MT-3 were installed at the bottom of two thermokarst ponds. Probe site MT-1, which is located in an area at the initial stage of a developing thermokarst pond, is more likely a water-filled depression due to thaw subsidence rather than a genuine thermokarst pond. Temperature probe MT-1 lies within a 1 × 2 m elongated pond which was 0.32 m deep ( Figure SM4f (teco_a_1949819_sm2254.zip)), while temperature probe MT-3 was also in an 8 × 10 m elongated pond which was 0.55 m deep ( Figure SM4h (teco_a_1949819_sm2254.zip)).

Figure 3.

Climate variability at Umiujaq over the period from 1926 to 2019 based on mean annual air temperature (MAAT) data from Kuujjuarapik, Inukjuak, and Umiujaq in Nunavik (Québec). Climate data for Kuujjuarapik and Inukjuak are from Environment Canada (2020).

Surface temperature statistics for eight characteristic ground surface conditions, as well as air temperature, are given in Table 1, from which ground surface conditions favourable or detrimental for maintaining permafrost can be identified. For instance, since the areas of permafrost mounds which are covered by frost boils, mosses and lichens have mean annual surface temperatures (MASTs) significantly below 0°C (Table 1), these surface conditions are favourable for maintaining permafrost. In contrast, topographic depressions having thick snow cover in winter (often resulting from higher vegetation cover such as shrubs and conifers which help retain snow), as well as depressions forming thermokarst ponds, have high mean surface temperatures close to or above 0°C which promote permafrost degradation.

On the graphs of mean daily surface temperature (MDST) as a function of mean daily air temperature (MDAT) in Figures 4 and 5, simple linear regressions forced to pass through the origin (MDST = 0°C; MDAT = 0°C) are shown for MDATs above and below 0°C. These regressions provided the thawing and freezing slopes, respectively (m_thawing and m_freezing in Figures 4 and 5 and in Table 1), which represent surface/air temperature ratios. These slopes were then used in the numerical model to predict the corresponding ground surface temperatures for the given ground surface conditions and for any given past or future mean daily air temperature using the following equations:

For instance, on the top of permafrost mound 2, which is covered by mosses and lichens, the freezing slope was 0.655 (Figure 4c), reflecting a close coupling between the air and ground surface temperatures. The ground is thus efficiently cooled down in winter which was favourable for maintaining permafrost. However, for the area invaded by tall shrubs with a thick snow cover in winter, the freezing slope was 0.0114 (Figure 4d), reflecting surface temperatures which stayed remarkably stable at 0°C, independent of air temperature. Field observations at the study site showed that permafrost does not form under such conditions. Moreover, for the two temperature probes MT-1 and MT-3 located in thermokarst ponds on the top of permafrost mound 2, the freezing slopes were 0.227 and 0.0401, respectively (Figure 5c and 5d). These thermokarst ponds froze to the bottom at some point during the winter. In the deep pond (probe MT-3), complete freezing occurred later in winter (in February 2019) compared to December 2018 for the shallow pond (probe MT-1). The freezing slope of temperature probe MT-1 was therefore greater than the slope of temperature probe MT-2. Before complete freezing, the bottom temperature was near 0°C due to the release of the latent heat of freezing of water. After freezing, the ground surface temperature cooled down due to heat extraction by conduction through the ice cover (Figure 5).

Figure 4.

Time series of hourly air and ground temperatures measured at selected probe sites: (a) TR-6 covered by mosses and lichens, and (b) TR-1 invaded by dense shrub vegetation, as a function of time from July 2018 to July 2019. Mean daily surface temperature (MDST) measured at the same probe sites (c) TR-6 and (d) TR-1 as a function of mean daily air temperature (MDAT) for the period from July 2018 to July 2019. Simple linear regressions through the origin (0, 0), used in the numerical modelling, are also provided. Location and photographs of temperature probe sites TR-1 and TR-6 are shown in Figures 1,  SM3 and SM4 (teco_a_1949819_sm2254.zip).

Except for the surfaces covered by vegetation, the thawing slopes for all other ground surface conditions, including the thermokarst ponds, were above 0.57, again reflecting a close coupling between the air and ground surface temperatures. Under a vegetation cover, the thawing slopes could be as low as 0.088 for the spruce and conifer coverage ( Figure SM8 (teco_a_1949819_sm2254.zip)), reflecting surface temperatures which stayed remarkably stable, being slightly above 0°C.

The time series of air and surface temperatures and corresponding graphs of MDST as a function of MDAT for the six other ground surface conditions listed in Table 1 are provided as supplementary online material ( Figures SM5 to SM10 (teco_a_1949819_sm2254.zip)).

Characterization of surface conditions

The field observations described above, including the infrared thermal imaging and surface temperature monitoring, allowed to define the following sequence from the most favourable surface condition to the most detrimental condition with respect to maintaining permafrost (Figure 1 and Table 1): (1) mosses and lichens, (2) frost boils, (3) development of shrubs and spruces, and (4) topographic depressions that eventually fill with water and gradually become thermokarst ponds (5). Characterizing these conditions was needed to choose the most relevant conditions to vary in the numerical model.

Figure 5.

Time series of hourly air and ground surface temperatures measured at selected probe sites in two thermokarst ponds: (a) the initial stage of thermokarst pond MT-1, and (b) the final stage of thermokarst pond MT-3, as a function of time from July 2018 to July 2019. Mean daily surface temperature (MDST) as a function of mean daily air temperature (MDAT) for the period from July 2018 to July 2019, measured at the same probe sites (c) MT-1 and (d) MT-3. Simple linear regressions through the origin (0, 0), used in the numerical modelling, are also provided. Location and photographs of temperature probe sites MT-1 and MT-3 are shown in  Figures SM3 and SM4 (teco_a_1949819_sm2254.zip).

As observed in the first available aerial photographs taken in 1957 (Fortier and Aubé-Maurice 2008), frost boils completely covered all permafrost mounds in the area of Umiujaq, being indicative of minimal snow cover and well-developed permafrost mounds. Due to the trend of climate warming recently observed in Nunavik from 1993 to 2003 (Figure 3), freezing action in the frost boils became less pronounced which allowed mosses and lichens to progressively colonize the top of these mounds. This subtle change in surface conditions was enough to form a thin insulating snow cover, which slightly increases the surface temperature during the winter (Table 1). However, the change in surface reflectivity in summer due to the mosses and lichens tends to decrease the summer surface temperatures (Figure 4a vs  Figure SM5a (teco_a_1949819_sm2254.zip)). Overall, the MAST of a surface colonized by mosses and lichens is lower than for frost boils which are free of vegetation and snow.

Thawing of ice-rich permafrost due to climate warming is accompanied by thaw consolidation and surface thaw settlement leaving localized depressions mainly located on the sides of permafrost mounds. Snow and water can accumulate in these depressions which promotes thermal insulation in winter and a latent heat effect of water, respectively, which increases the surface temperature, and induces more permafrost degradation. At some point, the localized depressions coalesce, fill with water, and a thermokarst pond forms. Finally, due to climate warming and related greening of the north (Myers-Smith and Hik 2017), including the development of shrubs and spruces on the less wind-exposed sides of the permafrost mounds, a thick insulating snow cover can accumulate over the vegetation which further contributes to permafrost degradation.

Table 1.

Statistics on ground surface temperatures for different surface condition characteristics of the study site and air temperature in the Tasiapik Valley at Umiujaq from 4 July 2018 to 3 July 2019 (MinST: minimum surface temperature; MaxST: maximum surface temperature; MAST: mean annual surface temperature; MST ≤ 0°C: mean surface temperature below or at 0°C; MST > 0°C: mean surface temperature above 0°C; MinAT: minimum air temperature; MaxAT: maximum air temperature; MAAT: mean annual air temperature; m_freezing: freezing slope; m_thawing: thawing slope). Surface conditions are listed in order from those more favourable to those more detrimental for maintaining permafrost.

The above-described characterization of ground surface conditions based on field observations of permafrost degradation may vary from one location to another. Since these changes in ground surface conditions are interrelated, a slight disturbance in the thermal regime of permafrost due to climate warming results in complex feedback loops where changes in ground surface conditions trigger permafrost degradation, which in turn affects the ground surface conditions. These changes in ground surface conditions, simulated as changes in surface temperatures, were used in the numerical simulations presented below to predict their impact on permafrost degradation. However, feedback from the simulated permafrost degradation back to the ground surface conditions was not considered.

Changes in three specific ground surface conditions were investigated using the numerical model: (1) a surface colonized by mosses and lichens which is the most favourable condition for maintaining permafrost, (2) shrubification, and (3) formation of a thermokarst pond which promotes permafrost degradation.

Numerical model development

Theoretical approach

In this study, a 2D numerical site model was developed using the HEATFLOW numerical code (Molson and Frind 2020) to better understand permafrost dynamics and response of permafrost to changes in ground surface conditions. The HEATFLOW code couples groundwater flow and advective-conductive heat transfer, including heat exchange at the air–ground interface, water-ice phase change, and latent heat within a heterogeneous and anisotropic environment. Fully saturated conditions are assumed for all simulations.

The general advective-conductive heat transfer equation which is solved in HEATFLOW is expressed as follows:

where θ is the porosity, S_w is the unfrozen water saturation, c_w is the specific heat of water (J/kg/K), ρ_w is the water density (kg/m³), v_i is the mean linear groundwater flow velocity (m/s), T is the temperature of the porous medium (°C), λ_B is the apparent thermal conductivity of the bulk porous medium (W/m/K) which is a function of λ_ice, λ_water and λ_solids and the respective water saturation, D_i,j is the hydrodynamic dispersion coefficient (m²/s), C_o is the effective heat capacity of the porous medium (J/ m³/K), x_i,j are spatial coordinates (m) and t is time (s). Flow velocities in Equations 3 are computed from a coupled transient density-dependent flow equation based on the conservation of fluid mass and Darcy's Law (Molson and Frind 2020).

In Equations 3, the effective heat capacity of the porous medium C_o is defined as:

where L is the latent heat of water (J/kg) and the subscripts w, i, and s refer to water, ice, and solids, respectively.

The unfrozen water saturation function W_u and the relative permeability k_r are expressed as follows:

where p is the terminal fraction of unfrozen moisture at a very low temperature, and q and Ω are empirical constants.

HEATFLOW uses the Galerkin finite element approach using rectangular prismatic elements with Picard iteration to solve the coupled groundwater flow and heat transport equations. Verification examples for various freeze-thaw problems are provided in Grenier et al. (2018) and Molson and Frind (2020). Field-scale applications to permafrost-impacted systems include those presented by Shojae-Ghias et al. (2018) and Dagenais et al. (2020).

Model domain and physical properties

Based on the conceptual model shown in Figure 1, a 2D vertical-plane numerical model was developed for permafrost mound 2 (see Figure 2c) which is the most instrumented of the two permafrost mounds at the study site. The model domain is 150 m long and 50 m in height (Figure 6), extending from the water table to the impermeable bedrock. The assumption of fully saturated conditions is considered is considered reasonable at the scale of 150 m and where the water table is within a couple of metres of ground surface.

The domain is divided into four distinct stratigraphic layers reflecting the local stratigraphy (Fortier et al. 2017, 2020) which is defined from top to bottom as: (1) a shallow sand aquifer, (2) a frost-susceptible silt aquitard, (3) a fluvio-glacial sand and gravel confined aquifer, and (4) impermeable quartzite bedrock (Figure 6). The initial body of ice-rich permafrost in the frozen silt layer extends horizontally from 60 to 100 m and from 22 to 45 m in elevation above sea level (elevation asl). Finally, the model is composed of a uniform mesh of 151 × 51 (= 7701) nodes in the horizontal and vertical directions, respectively, with all elements of dimensions 1 m × 1 m (see the grid on the right side of Figure 6).

Figure 6.

Numerical model domain based on the conceptual model of permafrost mound 2 (Figure 1), including the mesh (representative section shown to the right) and thermo-hydraulic properties of each layer (λ_s is the thermal conductivity of solids, K is the hydraulic conductivity, and θ is the porosity). The blue dashed line rectangle identifies the initial condition assumed for the permafrost block.

Hydraulic conductivities of all units were initially based on field data (Lemieux et al. 2020) and on a calibrated larger-scale numerical model of the study site (Dagenais et al. 2020). Conductivities were then further refined during calibration of the current smaller-scale model to observed temperature profiles. For the two principal units of marine silt and coarse sand and gravel, the final values (shown in Figure 6) were 10 times higher than those calibrated by Dagenais et al. (2020), and respectively 30 times lower and 15 times higher than reported by Lemieux et al. (2020); the latter, however, were based on only a few grain size analyses. All selected hydraulic conductivities are assumed isotropic. Solid phase thermal conductivities (λ_s) and porosities (θ) of all sediments were obtained from Dagenais et al. (2020) which are based on lab data and model calibration. For the quartzite bedrock, the assumed mean thermal conductivity of the solid matrix (λ_s) is 4 W/m/K (Cermak et al. 1982; Shim and Park 2013; Lee and Park 2015; Bédard et al. 2016) and the porosity (θ) is 0.038 (Hirth and Tullis 1989). The unfrozen moisture saturation and relative permeability functions for the silt layer, and their defining parameters from Equations 5a,5b,5c and 6, are provided in Figure 7.

Boundary and initial conditions

Flow system. Boundary conditions for the groundwater flow and heat transport systems are shown in Figure 8a and 8b, respectively. For the groundwater flow system (Figure 8a), the water table was represented as a quarter-cosine function with hydraulic heads decreasing to the right from 53 m to 50 m, consistent with the observed shallow water table. A hydraulic gradient was thus induced towards the discharge zone, reaching a maximum at the top-right model boundary representing the central creek in the Tasiapik Valley (Figures 1 and 2).

Figure 7.

Unfrozen water saturation W_u and relative permeability k_r of the silt layer as a function of temperature. The porosity θ of the silt is 0.4. Other parameters are defined in Equation 5 and 6

As a Type-1 boundary condition, the imposed water table allowed recharge to naturally vary depending on its frozen state and relative permeability. These conditions produced a reasonable maximum recharge in summer of the order 100 mm/yr which is consistent with that derived by Lemieux et al. (2020). The chosen water table gradient and quarter-cosine function were also confirmed to produce a Darcy flux in the confined aquifer on the order of 3.5 × 10^–7 m/s which is consistent with that simulated by Dagenais et al. (2020) and with a flux of 10^–6 m/s observed in a tracer experiment by Jamin et al. (2020).

Figure 8.

(a) Boundary conditions for the groundwater flow system of the numerical model of permafrost mound 2. Dominant flow directions in the shallow and confined aquifers are identified by black arrows. (b) Boundary and initial temperature conditions (T₀) for the corresponding heat transfer system. See Figure 6 for color coding of the units of the conceptual model.

At the left boundary, a hydraulic head of 53 m was imposed in the two aquifers, which corresponds to the water table elevation. All other flow boundaries at the right, left, and base of the model are assumed to be symmetric or impermeable.

An initial head of h_o = 53 m was assumed everywhere which rapidly equilibrated to the imposed initial and boundary conditions within the first few time steps of the spin-up period.

Heat transport system. Boundary conditions for the heat transfer system are shown in Figure 8b. Imposed (Type-1) temperatures were assigned along most of the upper boundary of the model which depend on the MAAT time series air temperatures at Umiujaq (Figure 3) and on the thawing and freezing indices (m_thawing and m_freezing) for different ground surface conditions as defined in Equations 1 and 2 (see also Figures 4 and 5, and Table 1). These imposed surface temperatures also depend on the local ground surface conditions which change in time.

A 10-m section at the far-right surface discharge zone was left as a zero-gradient temperature condition which allowed water to exit at its ‘natural’ transient temperature determined by the physical processes simulated in the model.

All remaining boundaries at the right and left faces of the thermal domain were assigned zero-temperature gradient conditions which allow advective heat transfer on the open parts of the left boundary and no thermal conduction along the corresponding impermeable boundaries. Finally, a geothermal heat flux of 0.032 W/m², characteristic of the Umiujaq region (Mareschal and Jaupart 2004), was imposed at the base of the model.

Time steps in the numerical model were kept at 1 day throughout the simulation time for each scenario. The head and velocity solutions were updated at each time step to account for temperature-dependent properties and to account for the changing air temperatures and ground surface conditions over time. Convergence criteria were 0.01 m for hydraulic head and 0.05°C for temperature which were usually reached within 2–3 iterations at each time step.

Figure 9.

Variations in mean annual air temperature from 1825 to 2100 used for the simulations: stable climate from 1825 to 1925, recent climate variability from 1926 to 2019 (see Figure 3), and assumed climate scenarios T1 and T2 from 2020 to 2100 (see Table 2).

Model calibration

Transient steady state in 2019

Simulated temperature distributions within the 2D section at the end of the spin-up period in January and August 2019 are shown in Figure 10a and 10b, respectively. The coldest conditions are found in January while the warmest conditions are found close to the surface in August when the active layer reaches its maximum depth. Simulated permafrost temperatures range from 0°C to –2°C, while ground temperatures surrounding the permafrost mound range from 2°C to 4°C.

Simulated and observed ground temperature profiles upstream, downstream, and within the permafrost mound are shown in Figure 10c and 10d, respectively. The corresponding simulated and observed profiles differ somewhat, with the simulated winter (January) temperatures generally warmer than those observed while the simulated summer (August) temperatures are lower than those observed. In all profiles, the simulated depth of zero annual amplitude is also somewhat greater than what was actually observed. Nevertheless, the general trends are suitably reproduced. The upstream temperature profile is warmer than the downstream profile (Figure 10c) due to the cold permafrost mound in between, which cools the groundwater as it flows downstream. Although the simulated ground temperatures are lower in 2019 than the available observed temperatures in 2018, the dimensions of the simulated permafrost mound relative to those of real permafrost mounds in the Tasiapik Valley are reasonable (Figure 10a and 10b). While the simulated thickness of the active layer of 1.8 m is slightly lower than the thickness of 1.9 m measured at the VDT-SYBU meteorological station (Figure 10b), the permafrost base at a depth of 22 m is similar in both simulated and measured cases. Simulated ground temperatures in the permafrost mound seem closer to those observed in 2001, where permafrost degradation was less significant and ground temperatures were colder than they are today. These differences in the thermal profiles are highlighted in the Discussion section but remain difficult to explain.

Figure 10.

Simulated ground temperatures and groundwater flow paths within the 2D domain at the end of the spin-up period in (a) January and (b) August 2019. (c) Simulated ground temperature profiles upstream and downstream, and the observed profile downstream of permafrost mound 2, in August and January 2019. (d) Simulated ground temperature profiles in permafrost mound 2 in August and January 2019, and measured profiles in August and January for both 2001 and 2018.

Scenarios of climate warming and changes in ground surface conditions from 2020 to 2100

Two scenarios of climate warming were considered in the predictive simulations from 2020 to 2100 (Figure 9): (1) scenario T1 extends the observed trend over the last 30 years with a linear increase in MAAT of 0.23°C from 2020 to 2100 or 0.003°C/yr (green dashed line in Figure 10), while (2) scenario T2 assumes an increase of 2°C over the same period or 0.025°C/yr (red dashed line in Figure 9) (Ouranos 2021).

Three scenarios of changes in ground surface conditions were also considered (Figure 11). In each case, the air temperature (Figure 9) is converted into imposed ground surface temperatures at the top boundary of the model based on the local ground surface condition and the corresponding thawing and freezing slopes in Equations 1 and 2 (Figures 4 and 5, and Table 1). Scenario SC1 is a reference scenario with no surface modification from 1926 to 2100 (Figure 11a). A case of increasing ground temperatures associated with the development of a thermokarst pond on the top of the permafrost mound is simulated in scenario SC2 (Figure 11b). In this scenario, an initial stage of development of a 1-m wide thermokarst pond is simulated from 2020 to 2060, then its width is increased from 1 (dimension of MT-1) to 6 m (mean dimension between MT-2 and MT-3 was chosen since it is a sudden change of surface conditions in time) in 2060. For scenario SC3, a step-wise shrubification is applied on either side of the permafrost mound (Figure 11c) from an initial growth stage over the 2020–2050 period followed by a second growth stage over the 2050–2100 period. The shrubification is assumed to suddenly invade 6 m on either side of the permafrost mound in 2020 and an additional 3 m in 2050, for a total of 18 m of shrubification over the 2020–2100 period.

Figure 11.

Changes in ground surface conditions in the heat transfer system (Figure 8b) used in the simulation scenarios (Table 2 and Figure 9) over the period from 2020 to 2100: (a) no change in ground surface conditions (ground surface condition SC1), (b) development of a thermokarst pond (surface condition SC2; MT-1 from 2020 to 2060 and MT-3 from 2060 to 2100), and (c) shrubification (change in ground surface condition SC3).

The hypothetical scenarios of climate warming coupled to changes in ground surface conditions (scenarios S1 to S4) are summarized in Table 2.

Predictive simulations with climate change from 2020-2100

The results of the predictive simulations showing impacts of climate warming and changes in ground surface conditions on permafrost degradation are presented in Figure 12 as spatial distributions of ground temperature differences (ΔT) with respect to the reference case scenario S1 after 100 years in January and August 2100. The reference case scenario S1 includes a moderate climate warming rate (climate scenario T1) of 0.23°C from 2020 to 2100 with no change in ground surface conditions SC1 (Table 2 and Figures 9 and 11). Comparisons are shown between scenarios S2 and S1 (Figure 12a), between scenarios S3 and S1 (Figure 12b) and between scenarios S4 and S1 (Figure 12c). Moreover, a comparison between the reference case scenario S1 in 2100 and the simulated ground temperatures at the end of the spin-up period in 2019 is provided as supplementary online material to show the impacts of a slight climate warming on permafrost degradation without any change in ground surface conditions.

Table 2.

Description of the hypothetical simulation scenarios from 2020 to 2100 as a function of climate scenarios identified as T1 and T2 and changes in ground surface conditions identified as SC1, SC2, SC3, and SC4. See Figure 9 for graphical representation of climate scenarios T1 and T2.

Figure 12.

Differences in ground temperature ΔT (°C) in the 2D model between scenarios (Table 2) in January (left side) and August 2100 (right side): (a) between scenarios S2 and S1, (b) between scenarios S3 and S1, and (c) between scenarios S4 and S1. The 0°C isotherm in 2100 of scenario S1 is identified by the white dashed line, while that of scenarios S2, S3, and S4 is identified by the solid white line.

As shown in Figure 12a, which compares only the effect of a higher rate of climate warming but with no change in ground surface conditions (comparison between scenarios S2 and S1), climate scenario T2 with a 2°C increase from 2020 to 2100 induces a slight thinning and a narrowing of the permafrost mound of less than 2 m. However, under the same climate scenario T2, the development of a thermokarst pond (change in ground surface conditions SC2) induces significant permafrost degradation (Figure 12b). In this comparison between scenarios S3 and S1, average temperatures increase by more than 0.5°C within the permafrost mound and by as much as 1°C to 1.5°C immediately beneath the thermokarst pond. Moreover, the 0°C isotherm of scenario S3, which follows the permafrost base, is 2 m closer to the surface than in scenario S1 where the shift is more pronounced beneath the thermokarst pond. A talik is also present beneath the thermokarst pond between depths of 1 and 5 m. The impacts of scenario S3 are mainly felt beneath the thermokarst pond and at the permafrost base.

Shrubification (change in ground surface conditions SC3) coupled with a climate warming of 2°C from 2020 to 2100 (climate scenario T2) also causes significant permafrost degradation (Figure 12c). In this case, the imposed winter surface boundary temperatures are warmer following shrubification due to the thermal insulation effect of snow trapped in the tall shrubs. As a result, the permafrost extent decreased by about 7 m on either side of the permafrost mound (below the shrubs), from 2020 to 2100 and ground temperatures increased by more than 3°C below the vegetation cover. Furthermore, the permafrost base is around 20 m deep in 2100 for scenario S1, while it is only 16 m deep for scenario S4. The increase in permafrost temperatures in areas not affected by any overlying shrubification in scenario S4 ranges from 1°C to 1.5°C in August and 0.5°C to 1°C in January. Scenario S4 thus has the most significant impact on permafrost dynamics and degradation and, under these conditions, permafrost is expected to eventually disappear.

A slight climate warming of 0.23°C from 2020 to 2100 with no change in ground surface conditions induces a slight narrowing and a thinning of about 2 m of the permafrost mound ( Figure SM11 (teco_a_1949819_sm2254.zip)). The upward migration of the permafrost base is slightly more significant on the right side (downgradient with respect to groundwater flow) than on the left side (upgradient).

Applying transient ground surface temperatures to constrain the numerical model

The approach adopted in this study of applying spatially and temporally variable ground surface temperatures, derived from field-observed surface/air temperature ratios for different ground conditions, does not explicitly reproduce the local thermal processes associated with these conditions. It does, however, catch the impacts of these conditions on permafrost dynamics. For instance, shrubification tends to trap more snow, thus increasing thermal insulation of the ground against cold winter air temperatures. In summer, shrubification increases the surface albedo, thus keeping the ground relatively cooler than bare soil. The impacts of these thermal processes on the ground surface temperature relative to the air temperature are visible in Figure 4b and 4d (see also  Figure SM8 (teco_a_1949819_sm2254.zip) for black spruce cover). In addition, a thermokarst pond can also favor snow accumulation and thermal insulation as with shrubification. Furthermore, the water column and ice cover in a thermokarst pond will delay the onset of winter freezing and spring thaw, respectively, due to the latent heat effect and thermal inertia (Figure 5). If the ice cover in a thermokarst pond is bottom-fast at some point in winter (ex. in a shallow pond), more rapid thermal conduction and ground cooling will occur (Figure 5) since the thermal conductivity of ice is significantly higher than water.

Since the numerical model developed herein also simulates heat advection, heat transfer due to groundwater flow at depth induces differences in ground temperature upstream and downstream of the permafrost mound (Figure 10). The local groundwater state and ground temperatures are also affected by changing recharge and discharge conditions, and by flow paths around the thermokarst pond (Figure 12b).

The simulation of these complex spatially and temporally-variable physical processes at the ground surface is challenging in field-scale numerical modelling and requires field measurement and monitoring of many parameters at high cost and effort. Nevertheless, since the applied ground surface temperatures used in the simulations presented herein already include the effect of these local-scale physical processes, the model is considered sufficiently constrained to accurately predict their impacts on permafrost dynamics while significantly reducing the model complexity and computational effort.

Comparison with other studies

A variety of approaches have been applied in cryohydrogeological modelling to simulate ground surface conditions, including the conceptual heat exchange layer and energy balance approaches (Lamontagne-Hallé et al. 2020). The results of this current study, using independent observed surface/air temperature ratios, provide similar conclusions but have provided a better understanding of the role of specific surface conditions on permafrost degradation. These results are consistent with those of Fisher et al. (2016) who identified important controls of ground surface conditions, in particular vegetation and soil characteristics (moisture), on the depth of the active layer, and with Albers et al. (2020) who used inverse modelling to show that simulated evolution of the thermal regime around a permafrost block was most sensitive to surface and near-surface ground conditions. McClymont et al. (2013) also concluded that differences in vegetation land cover have an important influence on thermal gradients. They showed that lateral thawing (3.5 m) of permafrost is faster than vertical thawing (1.4 m) over a period of 6 years (2004 to 2010) which is similar to our results, in particular for the shrubification scenario (Figure 12c; see section on Predictive simulations with climate change). Zhang et al. (2013) also demonstrated that it is crucial to consider changes in vegetation type in numerical modelling to better assess the impacts of climate warming on permafrost. Their model showed that the transformation of tundra to forest following permafrost thaw leads to a decreasing albedo which affects the energy balance of the lower atmosphere and tends to enhance the rate of climate warming and, therefore, the rate of permafrost thaw. Finally, Kurylyk et al. (2016) agree with the importance of considering land cover in numerical modelling and demonstrated that permafrost thaw enhances groundwater flow. Furthermore, with increasing precipitation from 10% to 25% until 2050 expected with climate change (Allard and Lemay 2012; Barrette et al. 2020), aquifer recharge will be greater, further enhancing permafrost degradation as shown by Douglas et al. (2020) and Mekonnen et al. (2021).

As discussed, permafrost evolution is driven by many local and regional factors that are interrelated. Numerical modelling remains the best way to anticipate physical processes related to permafrost degradation. Developing a complex numerical model considering all these factors would not necessarily generate a more useful model nor a better understanding (Batty and Torrens 2001; Hill 2006; Voss 2011; Brunetti et al. 2020), especially for sites without sufficient supporting data, as is common in northern groundwater environments.

In a continuous permafrost zone in north-eastern Siberia, Magnusson et al. (2020) showed that thermokarst ponds are characterized by a cyclic vegetation succession related to permafrost degradation and recovery. Younger thermokarst ponds are associated with dead shrubs which seems to increase thaw depths whereas older thermokarst ponds are associated with sedges and sphagnum which seems to show permafrost recovery. This is contrary to what the present study would have shown if the two surface conditions (i.e. shrubification and the development of a thermokarst pond) were combined in one model. Indeed, permafrost degradation in such a scenario would have been greater and the soil temperature warmer than with only shrubification or only development of a thermokarst pond. Differences between the simulations presented herein and the findings of Magnusson et al. (2020) could be due to the dynamics of snow cover as a function of vegetation cover, which would be difficult to simulate. Many consequences of climate change such as feedback loops between vegetation growth, snow cover dynamics, permafrost degradation, thaw settlement, and water-filled depressions, including those highlighted by Magnusson et al. (2020), are difficult to predict, even with advanced models. Although some responses to climate change are still unknown or not well understood, it is still expected that permafrost thaw is directly related to evolving ground surface conditions triggered by climate change and permafrost degradation.

Assumptions and limitations

Since models are never exact reproductions of reality, several inherent assumptions of the study presented herein need to be acknowledged. (1) Climate warming and changes in ground surface conditions are the most important factors controlling permafrost degradation. (2) Only two simple linear trends of climate warming were considered in the simulations while climate change will be much more complex in the future. (3) With climate change, precipitation is expected to increase in Nunavik (Zhang et al. 2019). However, since the imposed hydraulic gradient representing the water table was assumed to be fixed, changes in recharge were only related to freezing and thawing periods. (4) Applying transient ground surface temperatures is acknowledged as being an over-simplification of complex heat transfer processes that take place at the air/ground-surface interface. (5) The thawing and freezing slopes (m_thawing and m_freezing) for a given ground surface condition found from simple linear regressions on ground surface temperature as a function of air temperature are assumed constant over time. (6) Hysteresis effects in surface temperatures observed during ground freezing and thawing due to latent heat effects (Figure 4c, 4d, 5c and 5d) are not considered in the approach presented herein. This hysteresis effect could have been considered by producing daily ground-surface temperature vectors where, for the same air temperature, the surface temperature would depend on its timing relative to the transient freezing or thawing cycle. (7) Furthermore, the thawing and freezing slopes (m_freezing and m_thawing) were derived by forcing the linear regressions to pass through the origin, which would have somewhat altered their true slopes. Such forcing was needed to avoid any discontinuities in surface temperature when the air temperature crosses the 0°C isotherm either during cooling or warming. Nevertheless, neglecting hysteresis and applying average annual surface temperatures based on simple linear regressions as the surface boundary condition is considered a reasonable approach and was easy to implement in the numerical model for these long-term comparative simulations. (8) Numerical predictions of permafrost degradation were assumed sufficiently accurate under fully saturated conditions.

Several additional limitations in the modelling approach can also be highlighted. (1) The ground surface conditions of permafrost mounds at Umiujaq are highly variable, although for simplification, only the ground surface condition of mosses and lichens was considered during the spin-up period, while variations in surface topography were assumed negligible. (2) Changes in ground surface conditions observed in the field are gradual over time, while the model assumed abrupt changes at specific times in the future. (3) In simulating the impacts of changes in ground surface conditions, only one change at a time was considered whereas in the field these changes would occur concurrently. This strategy of changing only one surface condition at a time allowed a better understanding of the direct influence of each change on permafrost degradation. (4) Only two types of changes in ground surface conditions, shrubification and development of a thermokarst pond, were considered in the simulations presented herein while field conditions are much more complex with several types of changes occurring both spatially and temporally. Nevertheless, the simulation results presented herein still represent more complexity than what has usually been included in similar studies (see the previous section Comparison with other studies). (5) The model only considers a single homogeneous permafrost mound, whereas in the field, permafrost mounds would be heterogeneous, and those upstream of the study site would likely contribute to cooling of the groundwater temperature arriving at a downstream permafrost mound. At another site of discontinuous permafrost north of Umiujaq, Vallée and Payette (2007) show that the number and specific locations of permafrost mounds with respect to surface water can indeed affect degradation patterns. (6) With the exception of the ground surface boundary, the model boundary conditions were assumed constant in time, while in the field they could change. (7) Although the bedrock was assumed impermeable, fractures have been observed in this unit which could form preferential groundwater flow paths. (8) Feedback from permafrost degradation back to ground surface conditions is not considered. Rather, the imposed surface temperatures in the predictive scenarios are assumed reasonable estimates of the actual changes which would occur in the field. (9) Coarse vertical discretization with 1-m thick elements close to the surface may affect the simulation of heat transfer in the surficial layer where large variations in temperature both in time and depth occur. Using thinner elements in the superficial layer down to the depth of zero annual amplitude might have better reproduced the heat transfer processes close to the surface with a possibly better match between simulated and observed ground temperatures.

The above simplifications in the modelling approach could be justified since the objective of this study was to show the importance of considering changes in surface conditions over time in numerical modelling of permafrost dynamics and not to accurately reproduce the field observations. Nevertheless, the simulated extent and temperatures of permafrost mound 2 in 2019 are considered a reasonable match to the observed data. Including more detailed system geometry with a finer discretization at the surface of the model, combined with a full energy balance condition at the ground surface and additional calibration could help provide a better match between the simulation results and field observations. However, these improvements would not be expected to change the general conclusions regarding the relative importance of changes in ground surface conditions on permafrost degradation.