Numerical simulations were carried out based on a conceptual cryohydrogeological model of a permafrost mound near Umiujaq, Nunavik (Québec), Canada, to assess the impacts of climate warming and changes in surface conditions on permafrost degradation. The 2D model includes groundwater flow, advective-conductive heat transport, phase change and latent heat. Changes in surface conditions which are characteristic of the site were represented empirically in the model by applying spatially- and temporally-variable ground surface temperatures derived from linear regressions between monitored surface and air temperatures. After reaching a transient steady-state condition close to present-day conditions, the simulations were then extended to 2100 under hypothetical climate warming scenarios and using imposed changes in surface conditions consistent with observed on-site evolution. The simulations show that the development of a thermokarst pond and shrubification respectively induce ground warming of up to 0.5°C and 1.5°C, upward migration of the permafrost base by up to 2 and 4 m, and a decrease in the lateral permafrost extent of 1 and 7 m, relative to a reference case without changes in surface conditions. Feedback from permafrost degradation which drives changes in ground surface conditions should be included in future numerical modelling of permafrost dynamics.
Permafrost degradation is currently occurring in arctic and subarctic regions as a result of climate warming (Payette et al. 2004; Fortier and Aubé-Maurice 2008; Thibault and Payette 2009; Box et al. 2019). Evidence of permafrost degradation includes thickening of the active layer at the expense of permafrost, thaw settlement of ice-rich permafrost, formation of thermokarst lakes, appearance of taliks, and landslides occurring along permafrost slopes (Zhang et al. 2006; Fortier et al. 2008, 2011; Gaanderse et al. 2018). Moreover, as part of a feedback loop induced by global warming, the greening or shrubification of the north (Myers-Smith and Hik 2017) and local changes in topography due to thaw settlement promotes snow accumulation and thermal insulation which further contributes to permafrost degradation (Goodrich 1982; Roche and Allard 1996; Fortier et al. 2011). Fisher et al. (2016) specifically highlight the importance of ground surface conditions on the shallow thermal regime including the thickness of the active layer.
Permafrost thaw has significant impacts on the dynamics of northern ecosystems and on the service life of civil infrastructure in northern communities (Fortier et al. 2011; Pelletier et al. 2018). Understanding the complex physical processes leading to permafrost thaw in the current context of climate warming is therefore critical for developing adaptation methods.
The thermal regime of permafrost is mainly controlled by the geothermal heat flux and by heat transfer at the air/ground interface (French 2007; Riseborough et al. 2008). While geothermal heat fluxes at depth can be assumed locally constant, ground surface heat fluxes which depend on air-ground temperature gradients can vary significantly at the local spatial scale and over time. At northern latitudes, differences between air temperature and ground surface temperatures have been observed to vary from 1°C to 6°C depending on the type of vegetation and snow cover conditions (Gold and Lachenbruch 1973; Roche and Allard 1996). In many areas at high northern latitudes, the terrain is currently transitioning (‘greening’) from tundric to boreal due to climate change (Gamache and Payette 2005; Beck et al. 2015; Ropars et al. 2015; Pelletier et al. 2018; Xiao-Jing et al. 2021). To various degrees, other physical processes and subsurface characteristics also control permafrost dynamics, including thermal insulation from snow cover (Roche and Allard 1996; Zhang 2005; Rushlow et al. 2020), thermal properties of soil layers (Farouki 1981), groundwater recharge (Young et al. 2020), and groundwater flow (Bense et al. 2012; Frampton et al. 2013; Kurylyk et al. 2014; Grenier et al. 2018; Dagenais et al. 2020; Rushlow et al. 2020). In addition, local land surface controls include the slope, albedo, surface roughness, and soil saturation.
A conceptual model of a characteristic permafrost environment in the discontinuous permafrost zone in Nunavik (Québec), Canada, is shown in Figure 1 which highlights the physical processes controlling permafrost dynamics and degradation. This conceptual model is based on field observations of permafrost mounds in the Tasiapik Valley which drains into Lake Tasiujaq (formerly Lake Guillaume-Delisle) near the Inuit community of Umiujaq in Nunavik (Buteau et al. 2004; Fortier and Aubé-Maurice 2008; Fortier et al. 2008, 2020; Pelletier et al. 2018). Formed during frost heave in a unit of frost-susceptible silty sediments around 2000 cal years BP (calibrated years before present), the permafrost mounds in this valley appear as raised periglacial landforms with dome-shaped surfaces known as lithalsas (Pissart 2002). Their tops are dotted with frost boils which are free of vegetation (Roche and Allard 1996; Fortier and Aubé-Maurice 2008; Pelletier et al. 2018). Under such surface conditions, strong winds leave the tops of the permafrost mounds essentially free of snow in winter. Groundwater flows through a shallow aquifer in the active layer and in a deeper coarse-grained aquifer confined by the unit of silty sediments. Advective heat transport from groundwater flow through these aquifers plays an important role on permafrost dynamics and limits the lateral extent and thickness of the local permafrost mounds (Buteau et al. 2004; Dagenais et al. 2020). Climate change induces additional processes in the permafrost environment at Umiujaq including permafrost thaw, thaw settlement in the form of localized depressions, snow and water accumulation in depressions, and growth of mosses and lichens. Eventually, shrubification attracts additional snow accumulation which promotes further permafrost degradation (Pelletier et al. 2018). Ultimately, thermokarst ponds, which are water-filled depressions with remnant ramparts, replace these permafrost mounds after complete permafrost thaw (Fortier and Aubé-Maurice 2008).
Although they clearly play a fundamental role in permafrost dynamics, such climate-warming induced changes in surface conditions which are partly related to permafrost degradation, have been only rarely considered in numerical modelling of permafrost dynamics. McClymont et al. (2013) developed a 2D thermal-conduction model to show how surface conditions affected thermal evolution of a permafrost mound in the Northwest Territories, Canada. Their model included the effects of a peat plateau and surrounding bogs and fens at the surface boundary, but did not simulate thermal advection or future climate change. Kurylyk et al. (2016) developed a 3D numerical model of the same permafrost mound including variable surface conditions. They also showed the importance of landscape evolution and lateral heat transfer on permafrost thaw.
One of the major challenges in numerical modelling of cryohydrogeological processes in permafrost environments is to define realistic conditions at the upper boundary for simulating heat transfer at the air/ground interface (Lamontagne-Hallé et al. 2020). Heat transfer by conduction, convection and radiation at this interface is complex and requires a significant amount of different types of data, which have to be carefully monitored in the field, to constrain numerical models (Riseborough et al. 2008).
In this study, two-dimensional (2D) numerical simulations were performed to provide insight into the role of climate warming and associated changes in ground surface conditions on permafrost thaw. A well-instrumented permafrost mound at Umiujaq is used as the basis for the conceptual model. Climate warming is here represented by increasing air temperatures combined with associated increases in ground surface temperatures which depend on characteristic surface conditions at the site. Evolving vegetation and formation of a thermokarst pond area also considered. Infrared thermal imaging and temperature data from thermal probes buried a few centimeters below ground are used to derive specific ground surface conditions and temperatures at various locations on and around the permafrost mound. An empirical method based on simple linear regressions between observed ground surface and air temperatures for these different surface conditions is then used to control thermal conditions at the air/ground interface. The primary objective of this study is to investigate the effects on permafrost dynamics of changes in ground surface conditions induced by climate change, and to evaluate the importance of including such conditions in cryohydrogeological modelling.
The study site is located in a 2 km2 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.
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 (MAATUMIUJAQ) from 1926 to 2000 was assessed by averaging the MAATs at Kuujjuarapik and Inukjuak (MAATKUUJJUARAPIK and MAATINUKJUAK). This climate variability at Umiujaq (Figure 3), which controls permafrost dynamics, is considered in the numerical simulations presented and discussed in the next sections.
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)).
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 (mthawing and mfreezing 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).
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.
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.
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; mfreezing: freezing slope; mthawing: 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
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, Sw is the unfrozen water saturation, cw is the specific heat of water (J/kg/K), ρw is the water density (kg/m3), vi 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, Di,j is the hydrodynamic dispersion coefficient (m2/s), Co is the effective heat capacity of the porous medium (J/ m3/K), xi,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 Co 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 Wu and the relative permeability kr 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).
The numerical model was first calibrated to roughly reproduce the observed temperature profiles in 2019 in permafrost mound 2, then applied to simulate hypothetical scenarios of future climate warming and changes in ground surface conditions. To ensure the initial conditions did not unrealistically constrain the calibration, they were set sufficiently back in time to 1825, producing a transient (year-over-year) steady-state condition by 2019, both with respect to groundwater flow and heat transfer.
A two-step equilibration (or spin-up) approach was adopted for the model calibration: (1) the numerical model was first subjected to an initial 100-year spin-up period from 1825 to 1925 where air temperatures and surface conditions were kept constant such that the numerical simulation converged towards an equilibrium thermo-hydraulic regime, followed by (2) a second spin-up period starting from this equilibrium regime, where the model was run from 1926 to the calibration year 2019 taking into account the climate variability as reconstructed at Umiujaq (Figure 3). From this calibrated condition in 2019, different hypothetical scenarios from 2020 to 2100 were then simulated in which the MAAT was assumed to linearly increase and where changes in ground surface conditions, represented as changes in surface temperatures, were applied over time. These simulated changes in ground surface conditions were namely the development of a thermokarst pond on top of the permafrost mound and vegetation growth on its sides.
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).
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).
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).
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 ho = 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 (mthawing and mfreezing) 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/m2, 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.
The spin-up period for model calibration was divided into two intervals (Figure 9): (1) a constant MAAT of –6.3°C with an amplitude of 19°C was applied from 1825 to 1925 assuming a cold stable temperature characteristic of the end of the Little Ice Age, and (2) using observed MAATs and the same annual amplitude over the period from 1926 to 2019 (Figure 3). At the start of the spin-up period, uniform initial temperatures were imposed for each unit (Figure 8b): 7°C for the surficial sand and silt layer, 10°C for the sand and gravel layer and in the deeper bedrock, and –2°C for the permafrost body.
For the spin-up period from 1825 to 1925, the ground surface conditions were assumed constant in time. Assigned surface temperatures were associated with mosses and lichens over the top 40-m of the permafrost mound and low shrubs on the flanks (Table 1 and Figure SM6 (teco_a_1949819_sm2254.zip)). In 1925, surface conditions on the flanks were changed to dense shrubs with a transition zone between the mosses and lichens and the low shrubs (Table 1 and Figure 4b and 4d) (Fortier and Aubé-Maurice 2008). The 40-m wide zone of mosses and lichens on the top of the permafrost mound did not change over the period from 1825 to 2019.
The simulated groundwater flow and temperature conditions at the end of the spin-up period in 2019 were used as the reference-year initial conditions for future predictions starting in 2020. Hypothetical scenarios of climate warming and changes in surface conditions were then simulated for the period from 2020 to 2100.
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.
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.
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.
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.
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 (mthawing and mfreezing) 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 (mfreezing and mthawing) 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.
Field investigations on and around permafrost mounds at Umiujaq, Nunavik (Quebec), Canada, support the following classification of impacts of changes in ground surface conditions on permafrost dynamics from those more favourable to those more detrimental for maintaining permafrost: (1) frost boils, and mosses and lichens on the top of permafrost mounds are most favourable, (2) localized thaw settlement allows snow cover accumulation and sporadic accumulation of surface water which tends to warm and thaw permafrost, (3) formation of thermokarst ponds when thaw settlement becomes significant is detrimental, and (4) shrubification of the sides and tops of permafrost mounds induces major and irreversible permafrost degradation.
The simulation results presented in this study show the importance of considering surface conditions in cryohydrogeological modelling and provide a better understanding of physical processes associated with permafrost degradation. Although the numerical modelling approach of applying transient ground surface temperatures as a boundary condition does not explicitly reproduce the local thermal processes associated with different ground surface conditions, the impacts of these conditions on permafrost dynamics and degradation can still be assessed. Comparing simulations with and without changes in ground surface conditions, the worst-case coupled scenario of climate warming by 2°C from 2020 to 2100 and a step-wise 18-m shrubification of a permafrost mound over the same period is characterized by ground warming up to 1.5°C, an upward migration of the permafrost base by up to 4 m, and a decrease in the extent of the permafrost mound of 7 m. Changes in ground surface conditions over time have a significant impact on heat transfer at the air-ground interface and disrupt the subsurface thermal regime. The numerical model developed in this study has led to a better understanding of the relationships between permafrost degradation and climate warming which induce changes in ground surface conditions, as represented by associated changes in surface temperatures. The influence of groundwater flow on permafrost dynamics was also shown.
Recommendations for future numerical modelling of permafrost dynamics and degradation include: (1) adding hysteresis effects due to phase change and latent heat in the relation between the air and ground-surface temperatures, (2) adding any additional upstream permafrost mounds to more accurately simulate the temperature of inflowing groundwater, (3) accounting for changes in precipitation over time in order to better simulate the impacts of climate warming, and (4) including direct feedback mechanisms between the degrading permafrost and changing ground surface conditions.
The authors wish to thank the Inuit community of Umiujaq for their kind hospitality and interest in this project. Logistical support from the Centre d'études nordiques (CEN) of Université Laval, over many years of field investigation and data collection, is gratefully acknowledged. We also thank Jonathan Fortin for his help in the field, and Guillaume Allard and Pierre Therrien of the Département de géologie et de genie géologique of Université Laval for technical help with computational issues. Thanks are also due to two anonymous reviewers and Hugo Asselin, the Editor-in-chief of Écoscience, for their constructive comments which improved the quality of the paper.
Financial support to the first author in the form of scholarships was provided through the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grants of the second and third authors. Field campaigns at Umiujaq were supported in part by a NSERC Strategic Grant, by funding from the Ministère de l'Environnement et de la Lutte contre les changements climatiques du Québec (MELCC), and by a grant from the Northern Scientific Training Program (NSTP) of Polar Knowledge Canada to the first author.