Articles | Volume 17, issue 8
Research article
24 Aug 2023
Research article |  | 24 Aug 2023

Investigating the thermal state of permafrost with Bayesian inverse modeling of heat transfer

Brian Groenke, Moritz Langer, Jan Nitzbon, Sebastian Westermann, Guillermo Gallego, and Julia Boike

Long-term measurements of permafrost temperatures do not provide a complete picture of the Arctic subsurface thermal regime. Regions with warmer permafrost often show little to no long-term change in ground temperature due to the uptake and release of latent heat during freezing and thawing. Thus, regions where the least warming is observed may also be the most vulnerable to permafrost degradation. Since direct measurements of ice and liquid water contents in the permafrost layer are not widely available, thermal modeling of the subsurface plays a crucial role in understanding how permafrost responds to changes in the local energy balance. In this work, we first analyze trends in observed air and permafrost temperatures at four sites within the continuous permafrost zone, where we find substantial variation in the apparent relationship between long-term changes in permafrost temperatures (0.02–0.16 K yr−1) and air temperature (0.09–0.11 K yr−1). We then apply recently developed Bayesian inversion methods to link observed changes in borehole temperatures to unobserved changes in latent heat and active layer thickness using a transient model of heat conduction with phase change. Our results suggest that the degree to which recent warming trends correlate with permafrost thaw depends strongly on both soil freezing characteristics and historical climatology. At the warmest site, a 9 m borehole near Ny-Ålesund, Svalbard, modeled active layer thickness increases by an average of 13 ± 1 cm K−1 rise in mean annual ground temperature. In stark contrast, modeled rates of thaw at one of the colder sites, a borehole on Samoylov Island in the Lena River delta, appear far less sensitive to temperature change, with a negligible effect of 1 ± 1 cm K−1. Although our study is limited to just four sites, the results urge caution in the interpretation and comparison of warming trends in Arctic boreholes, indicating significant uncertainty in their implications for the current and future thermal state of permafrost.

1 Introduction

It is well known that the climate is changing, and the impact of climate warming on polar regions is disproportionately severe due to Arctic amplification (Serreze and Francis2006). Arctic regions are warming 2–3 times faster than the global average (Isaksen et al.2022; Hinzman et al.2005) which in turn is causing accelerated degradation of permafrost (Jorgenson et al.2006). Degrading permafrost can have a substantial impact on local hydrology and ecological systems (Schuur and Mack2018) through the formation of thermokarst lakes, the formation of trough networks, and landscape deformation. In addition, rapid changes in landscape can also pose a threat to human populations in Arctic regions by destabilizing existing infrastructure (Schneider von Deimling et al.2021; Nelson et al.2002). Understanding the past, present, and future evolution of Arctic permafrost therefore plays a crucial role in assessing, predicting, and mitigating the impacts of climate change.

A significant challenge in monitoring permafrost thaw, however, is the general lack of high quality and accessible long-term observational data in the Arctic. Over the past few decades, circumpolar borehole temperature measurement sites have been established to monitor changes in permafrost environments. Many of these sites include borehole sensor arrays which provide measurements of ground temperature as deep as 50 m or more below the surface. Recent studies have attempted to leverage these borehole data from the Global Terrestrial Network for Permafrost (GTN-P), a database of shallow and deep ground temperature measurements taken from borehole sites across the Arctic (Biskaborn et al.2015; Burgess et al.2000), to quantify changes in permafrost temperatures at a global scale (Chen et al.2021; Biskaborn et al.2019; Romanovsky et al.2007). These studies have consistently found that permafrost is warming globally (0.01 to 0.06 K yr−1 in the continuous permafrost zone) in response to the rapidly changing Arctic climate.

Temperature, however, describes only one part of the thermal state of permafrost (Smith et al.2022). A considerable fraction of energy can be consumed and released during thawing and freezing in the form of latent heat, thereby masking trends in the temperature signal within the critical range where phase change occurs (Nicolsky and Romanovsky2018; Romanovsky et al.2010; Romanovsky and Osterkamp2000; Riseborough1990). This partitioning of energy can be expressed mathematically as the sum of two constituent energy contents:

(1) H = H S + H L ,

where HS=TC(θ(T)) and HL=Lθ(T) [J m−3] are referred to as sensible and latent heat, respectively. Here T [K] is temperature, θ(T) [–] is (temperature dependent) volumetric liquid water content, C(θ) [JK-1m-3] is volumetric heat capacity, and L [J m−3] is the latent heat of fusion of water. The relationship between sensible and latent heat, and thus consequently the rate of permafrost thaw, depends heavily on the nonlinear freezing characteristic θ(T), often referred to as the soil freezing characteristic curve (SFCC) (Koopmans and Miller1966) for porous material. The SFCC governs the nonlinear relationship between liquid water content and temperature during phase change, which arises due to the capillary pressure difference between frozen and unfrozen water present in the pore space. The functional form of the SFCC is empirically determined from soil properties, which typically requires measurements of soil moisture content over the critical temperature range. However, such measurements are not available at most borehole sites. The resulting lack of information on the distribution of energy in the ground represents a significant source of uncertainty in understanding what temperature trends observed in many GTN-P boreholes really tell us about the current thermal state of permafrost across regions. In order to understand how permafrost is changing, it is vital that we find ways to quantify this uncertainty and develop physically coherent models to accurately link observed temperature trends to unobserved or sparsely observed quantities such as latent heat and active layer thickness (ALT).

The problem of linking unobserved components of the thermal state to observed quantities (e.g., temperature) is fundamentally a problem of inversion or inference under uncertainty. In other words, we wish to recover the full thermal state of permafrost (including the distribution of latent heat) given observed borehole temperatures and some partial, a priori information about the local landscape and meteorological conditions. Bayesian inference solves this problem probabilistically by combining information from observations with prior knowledge in order to obtain a distribution over the unobserved variables. These methods are particularly attractive for geoscientific applications where observational data often come from systems that are neither stationary nor well controlled (e.g., the local environment or the broader Earth system), thereby making traditional, frequency-based statistical inference difficult to apply and interpret correctly (Caers2018). As a result, Bayesian analyses have become increasingly common in the geosciences (Caers2018; Qu et al.2008; Berliner2003; Franks and Beven1997) and, more recently, also in cryospheric research (Gregory et al.2021; Garnello et al.2021; Verjans et al.2020; Gopalan et al.2018; Wainwright et al.2017).

Unfortunately, fully-fledged Bayesian inference using standard numerical sampling methods is often infeasible for complex, simulation-based models where forward evaluation of the model over longer time periods is computationally expensive, which is the case for most transient models of dynamical processes (Cranmer et al.2020). A wide range of strategies for overcoming this challenge have been proposed. One common family of methods known as “approximate Bayesian computation” (ABC) attempt to circumvent the issue of computing the likelihood by means of rejection sampling and summary statistics (Sisson et al.2018; Rubin1984). Particle filter (also known as “sequential Monte Carlo”) methods (Sisson et al.2007; Doucet et al.2001) have also been widely employed due to their efficiency in solving nonlinear filtering problems (Kantas et al.2014; Noh et al.2011; Moradkhani et al.2005). Other more recently proposed methods include Bayesian evidential learning (BEL) (Thibaut et al.2022; Hermans et al.2018) and the “calibrate, emulate, sample” algorithm of Cleary et al. (2021), both of which leverage the predictive power of machine learning algorithms to obtain an approximate form of the posterior distribution using a small number of initial forward simulations and a learned mapping to the relevant predictors.

Figure 1Map of the four study sites overlaid on the permafrost extent from Obu et al. (2019). Orange marks the tree line (Brown et al.1997). Photos are reprinted/adapted with permission from published documentation for each site: Barrow (Bar2021), Samoylov Island (Boike et al.2019), Bayelva (Boike et al.2018a), and Parson's Lake (Smith et al.2018) (permission granted by Open Government License - Canada:, last access: 22 August 2023).

These computational challenges in inverse modeling are naturally also relevant for transient models of permafrost processes which typically require solving one or more partial differential equations governing heat and water flow in the forward evaluation (Westermann et al.2023; Tubini et al.2021; Riseborough et al.2008; Nicolsky et al.2007). Finding computationally feasible methods for performing simulation-based inference is, therefore, a key methodological challenge in investigating the thermal state of permafrost with numerical models. However, relatively few studies to date have attempted to address this challenge. Romanovsky and Osterkamp (1997) used an analytical equilibrium model to invert air and soil temperatures measured at three sites in northern Alaska; Nicolsky et al. (2009) later used a traditional variational approach to invert measured soil temperatures with a numerical model of 1D heat conduction similar to the one used in our study. Harp et al. (2016) analyzed the effects of soil property uncertainties on permafrost thaw using the null-space Monte Carlo (NSMC) method applied to the Advanced Terrestrial Simulator (ATS) of Coon et al. (2019). More recently, Garnello et al. (2021) used Bayesian methods to calibrate the GIPL permafrost model (Jafarov et al.2012) and make long-term permafrost thaw projections.

The primary objectives of this study are (i) to analyze the relationship between near-surface air and permafrost temperature trends observed at several research stations in the Arctic continuous permafrost zone and (ii) to relate observed borehole temperature trends to plausible changes in latent heat over the last 2 decades via ensemble inversion of a transient heat conduction model with phase change. In both cases, we focus on the probabilistic quantification of uncertainty using Bayesian inference. To overcome the aforementioned computational challenges in solving the inverse problem for heat conduction, we apply a recently developed particle-based sampling method for approximate posterior inference (Garbuno-Inigo et al.2020) suitable for complex numerical models. We apply our methods to four different sites (Fig. 1) located in different climate zones in northeastern Siberia, Svalbard, northern Alaska, and northwestern Canada, where long-term borehole temperature records (> 10 years) as well as stratigraphic information are available.

Boike et al. (2019)Boike et al. (2018a)ECCC (2011)Theisen (2022)

Table 1Summary of climatological conditions for each of the four study sites included in this work. Air temperatures (for January and July) are monthly averages, rainfall refers to total annual unfrozen precipitation, and max snow depth refers to the maximum annual snow depth. Data are averaged over the reference period. All quantities were computed directly from published data for the indicated time period unless otherwise specified. Climatology data for Parson's Lake are taken from the nearest available long-term observation station at Trail Valley Creek, which is about 37 km southeast of the borehole site.

a  Computed using ECCC data from the later time period of 2015–2021. b Estimated visually from Fig. 5 of Burn and Kokelj (2009) for the earlier time period of 2000–2007. c Estimated visually from Fig. 6 of Shiklomanov et al. (2010) for the earlier time period of 1994–2010.

Download Print Version | Download XLSX

2 Study sites

In order to cover a broad range of environmental and climatic conditions, we apply our inverse modeling method to borehole and meteorological data obtained from sites in four different regions within the continuous permafrost zone. Figure 1 shows the geographical locations of each site along with photos of the local landscape. Table 1 provides a summary of available data on the long-term climatology of each site for a given reference period.

2.1 Samoylov

The Samoylov research station (Boike et al.2019) lies on Samoylov Island (722211.76′′ N, 1262851.0′′ E) within the Lena River delta at the edge of the Laptev Sea in northeastern Siberia. The region is characterized by ice-wedge polygonal tundra, dense vegetation coverage, large water bodies, and sandy ice-rich sediment. The site lies within the continuous permafrost zone and has permafrost reaching depths of 400 to 600 m below the surface (Grigoriev1960) with active layer thickness ranging from 0.3 to 0.7 m. Recent surface temperature reconstructions show that the permafrost has been subjected to substantial warming over the last century (Kneier et al.2018). The station provides automated monitoring of surface and climate conditions in addition to an instrumented 27 m borehole drilled in 2006 from which hourly ground temperature measurements are available up until September 2021.

2.2 Bayelva

The Bayelva station (Boike et al.2018a) is located on the Brøgger peninsula of the island of Spitsbergen, Svalbard, approximately 3 km from the small research village of Ny-Ålesund (785514.987′′ N, 11503.156′′ E). The West Spitsbergen Current helps to maintain comparatively mild seasonal temperatures in the area around Ny-Ålesund. The unglaciated coastal areas are known to have continuous permafrost extending to depths as deep as 100 m (Humlum2005), with active layer depths of 1 to 2 m. The site is situated on top of a hill, and the surrounding landscape is characterized by patterned permafrost ground and sparse (50 %–60 %) vegetation coverage. The soils on the hill range from silty loam to silty clay with some marine deposits, while coarser-grained silty sands can be found in the surrounding area. The site includes an instrumented 9 m deep borehole with hourly measurements of ground temperatures dating back to 2009.

2.3 Parson's Lake

The Parson's Lake, “T5 Upland”, borehole site (Smith et al.2018; Wolfe et al.2010) is situated east of the Mackenzie delta, northeast of Inuvik in the Northwest Territories and within the continuous permafrost zone of northern Canada (685729.879′′ N, 1335017.159′′ W). The region is classified as dwarf birch tundra with dense vegetation coverage and large water bodies (Burn and Kokelj2009). Soils in the area are mostly coarse-grained sand in the upper layers, with finer-grained silty-clay soils in the deeper layers (>1m). Soil cores also show evidence of excess ground ice in moderate amounts (approximately 5 %–25 % volumetric content) 1 to 2 m below the surface. Permafrost extends down to 100 m or more in the outer delta area (Allen et al.1988), and active layer thickness typically ranges from 0.35 to 0.65 m. The site features a roughly 10 m deep borehole with daily mean ground temperature measurements dating back to 2006.

2.4 Barrow

The Barrow North Meadow Lake 2 (NML2) borehole site (711914.07′′ N, 1563857.587′′ W) is situated near the village of Utqiaġvik (formerly known as Barrow), which lies at the northernmost point of Alaska along the Arctic coast. The surrounding landscape is characterized by drained thaw lakes and upland tundra consisting of gravel and silty soils with highly variable moisture and organic content in the active layer (Shiklomanov et al.2010). The region lies within the continuous permafrost zone and is known to have permafrost extending down to 300 m (Brewer1958) or more below the surface, with highly variable active layer thickness ranging from 0.1 to 1.0 m. The upper permafrost layer is fairly ice rich, with ground ice content as high as 50 %–75 % by volume in some areas. In this work, we use data from a borehole at the North Meadow Lake 2 (NML2) site, where instrumented daily average soil temperature measurements in the permafrost layer down to 15 m are sparsely available for the time frame 2012–2017, and annual measurements down to 47 m are available from 2009–2021 (Bar2021; Nelson et al.2008; Romanovsky et al.2002). We were only able to make use of the annual measurements below the assumed depth of zero annual amplitude in this study due to a large number of data gaps in the instrumented measurements.

3 Methods

3.1 Bayesian inference

The Bayesian approach to statistics provides a natural framework for inferring unobserved quantities of interest while simultaneously accounting for their associated uncertainties (Berliner2003). This is accomplished via Bayes' rule:

(2) p ( X | Y ) = p ( Y | X ) p ( X ) p ( Y )     with     p ( Y ) = x X p ( Y | X = x ) p ( X = x ) d x ,

which can be seen as a generic formula for obtaining the so-called posterior distribution of an unobserved quantity X given observations Y. The prior distribution p(X) reflects our uncertainty about X a priori (i.e., before observing Y), while the likelihood p(Y|X) measures how well the model's predictions agree with the observations Y. In this work, Y values are temperature measurements, typically sampled over time and/or space, whereas X values are unknown model parameters or unobserved physical quantities such as soil properties, active layer thickness, or the ratio of sensible to latent heat. The overall objective is then to obtain the posterior distribution p(X|Y) of these unknown parameters given the temperature measurements which quantifies not only the best-fitting parameter settings but also the associated modeling uncertainties.

Difficulties in the practical applications of Bayesian methods have historically arisen from the intractability of the integral in the denominator of Eq. (2), often referred to as the marginal likelihood p(Y). Advances over the last few decades in numerical sampling methods such as Markov Chain Monte Carlo (MCMC) (Hoffman and Gelman2014; Bishop2006; Duane and Kogut1986) which sidestep the need to compute the marginal likelihood, as well as a general increase in available computing power, have made Bayesian methods significantly more accessible. However, computing the likelihood still requires running the forward model given each set of parameters; for many scientific and simulation-based models, this incurs a significant computational burden.

3.2 Trend analysis

As a preliminary step of our study, we analyze trends in mean annual air and permafrost temperatures for each study site. For Samoylov, Bayelva, and Barrow, we use daily averaged in situ air temperature data where available and bias-corrected ERA5-Land (Mu noz-Sabater et al.2021) 2 m air temperature in all years in the study period (2000–2020) where in situ observations are unavailable. The bias-correction procedure for air temperature closely follows the empirical quantile mapping method of Piani et al. (2010), in which the empirical quantiles of both the model (reanalysis) and observational data are computed over some reference period (in this study, we use the full time period for which observations are available); the model data are then mapped to the corresponding quantiles of the observations. We deviate slightly from this procedure by only applying the correction to values within the observed (in situ) temperature range to avoid artificially clipping extreme temperatures in the reanalysis product. For the Parson's Lake borehole site, we use air temperatures from ERA5-Land without any correction since we were not able to find long-term in situ measurements near the station.

Due to the small sample size of our annual data (fewer than 20 years at all sites), we analyze air and permafrost temperature trends with a Bayesian robust linear trend model. We choose the Student's t distribution as the likelihood due to its longer tails which allow for more tolerance of unusually cold or warm years without severely affecting the overall trend (Lange et al.1989) . We use standard, weakly informative prior distributions for each parameter, with the exception of the slope μ1 to which we assign a mildly informative unit Gaussian. This is justified by the fact that annual average change in temperature near and below the surface can be reasonably expected to fall well below 1 K per annum. The specification of the full probability model can be found in Appendix B1.

Figure 2Illustrations of key concepts employed by our heat transfer model: (a) examples of soil freezing characteristic curves used in our model for two different soil textures. The curves shown in here are based on observations from Samoylov Island (sandy silt) and Bayelva (silty clay). (b) Simplified graphical representation of a four-layer stratigraphy configuration of the transient heat conduction model. The model uses stacked, homogeneous soil layers, where each layer is represented as a partitioning of excess ice, pore space, and solid material. The actual number of soil layers in the stratigraphy varies between sites. Tsurf and Qgeo represent the surface temperature and geothermal heat flux boundary conditions, respectively. Each ξ represents the characteristic fraction of one soil component.


3.3 Transient heat conduction model

3.3.1 Subsurface heat transfer

In order to simulate the diffusion of heat in permafrost soils, we use a transient heat conduction model based on the same numerical procedures and parameterizations of CryoGrid (Westermann et al.2023, 2016), which solves the 1-dimensional nonlinear heat equation with phase change:

(3) H ( T ) t = z k ( T ) T ( z , t ) z ,

where H(T) [J m−3] is the discrete enthalpy function over d finite volumes Δi(z), i=1,,d is defined on the depth domain z∈𝒵 with Z=[0m,1000m], and k(T) [Wm-1K-1] is the (temperature dependent) thermal conductivity. Note that Eq. (3) is in so-called mixed form, i.e., the conserved quantity (enthalpy) is linked to the diffused quantity (temperature) by a nonlinear constitutive relationship. We can recover temperature from enthalpy via the following nonlinear system:

(4) T = H - L θ ( T ) C ( θ ( T ) ) ,

which can be solved using standard numerical methods, such as Newton iteration, or approximated using a pre-calculated interpolant. Here, θ(T) represents the SFCC for a given soil type. In this work, we use the SFCC specified by Eqs. (21)–(23) in Dall'Amico et al. (2011) (see Fig. 2a and Appendix B2.1 for details), which is derived from the widely used van Genuchten (1980) soil-water retention equations.

Since no general solution exists for Eq. (3) with arbitrary boundary conditions, numerical approximation is necessary in order to find a solution for T(z,t). We discretize the spatial domain according to a fixed, non-uniform grid with the smallest grid cell spacing of 0.05 m from the surface down to a depth of approximately 4 m and a pseudo-exponentially increasing grid spacing thereafter.

We specify the soil stratigraphy for each site as stacked homogeneous layers, where the soil in each layer is parameterized according to the four characteristic fractions: excess ground ice (ξxic), natural porosity (ξpor), saturation (ξsat), and organic content (ξorg). From these independent fractions, the volumetric soil constituent fractions for excess ground ice (θx), pore water and ice (θwi), air (θa), and mineral and organic soil fraction (θm,θo) can be computed via recursive partitioning of the volume, as shown in Fig. 2 and described in Appendix B2.1. Note that, while excess ground ice is represented through super-saturation of the soil volume (effectively just increasing the porosity), soil and water displacement from ground ice melt and hydraulic processes are neglected, and the total water and ice content θwi+θx is assumed to remain constant throughout the simulation.

3.3.2 Initial and boundary conditions

Accurate simulation of the surface energy balance, snow dynamics, and surface hydrology is notoriously difficult and comes with many sources of uncertainty in landscapes with very heterogeneous surface characteristics (Langer et al.2011a, b). An empirical alternative to approximate the bulk thermal effect of these complex surface processes is to use n factors (Riseborough et al.2008; Lunardini1978), which attempt to approximate the long-term deviation of surface temperature from air temperature via lumped scaling factors:

(5) n F = FDD soil FDD air , n T = TDD soil TDD air , FDD = t 1 t 2 T - ( t ) d t , TDD = t 1 t 2 T + ( t ) d t ,

where T+(t)=max(T(t),0C), and T-(t)=min(T(t),0C). Note that FDD and TDD represent the commonly known freezing and thawing degree days. The resulting dimensionless factors nF,nT[0,1] can be used to approximate insulated surface temperatures by scaling the instantaneous air temperature:

(6) T surf ( t ) = n F ( t ) T air ( t ) T air ( t ) 0 n T ( t ) T air ( t ) T air ( t ) > 0 .

We use Tsurf, with Tair derived from the same, bias-corrected ERA5-Land air temperature record described in Sect. 3.2, as the Dirichlet upper boundary condition for temperature in the transient heat conduction model. In order to account for possible decadal scale changes in surface conditions, the seasonal n factors nF and nT are reparameterized as time-dependent, piecewise linear functions with knots at five pre-defined time points: 2000, 2005, 2010, 2015, and 2020. The values at each knot, e.g., nF(0), nF(1), nF(2), nF(3), and nF(4) (same for nT), are varied as unknown parameters, thereby allowing Tsurf, which is not generally known or measured reliably, to be adapted according to the observed borehole temperatures. The geothermal heat flux at the lower boundary (1 km below the surface) is set to 0.053 W m−2 and can generally be assumed to have a negligible effect on the decadal timescales considered in this work (Hermoso de Mendoza et al.2020).

The initial temperature profile at each site is parameterized similarly to the n factors as a piecewise linear function between five knots: 0m,T^air(t0)), (zALT,T1), (zZAA,T2), (zbase,0C), and (1000 m,T3), where T^air(t0) is the average annual air temperature for the first simulation year, and T1, T2, T3, and zbase (i.e., depth of the permafrost base) are varied as unknown parameters. zALT and zZAA nominally refer to the initial active layer thickness and depth of zero annual amplitude, respectively; these values are set based on rough estimates from the relevant literature for each site. Each simulation additionally includes a 20 year “spin-up” period (1979–1999) which is excluded from the analysis. The spin-up period allows the thermal state of the upper subsurface layers to reach dynamic equilibrium with the forcing, thereby minimizing the impact of the initial condition on the upper 10 m of the soil profile. Conversely, the parametric initial state prior to the spin-up period broadly accounts for first-order uncertainty in the initial temperature profile of the deeper layers which would otherwise require much longer spin-up periods (centuries or millennia) to equilibrate.

Figure 3ISO standard flowchart showing a simplified overview of the inversion workflow using ensemble Kalman sampling (EKS).


3.4 Inverse modeling

3.4.1 Theoretical framework

The inverse problem of recovering some unobserved sequence of states x1:τ from observations y1:τ over τ temporal samples can be expressed via Bayes' rule (Evensen et al.2022):

(7) p ( x 1 : τ , ϕ | y 1 : τ , s 1 : τ ) p ( y 1 : τ | x 1 : τ , ϕ , s 1 : τ ) × p ( x 1 : τ | ϕ , s 1 : τ ) p ( ϕ ) ,

where ϕ∈Φ represents the parameters characterizing the underlying hydrothermal processes in a dynamical model of the state variables x, and s1:τ represents any set of relevant exogenous variables (e.g., meteorological forcings, energy sources or sinks), which define the boundary conditions of the system. The objective of the inverse problem is then to obtain samples from the posterior distribution p(x1:τ,ϕ|y1:τ,s1:τ) over model states x1:τ and parameters ϕ. We assume that the parameters ϕ are independent of the exogenous inputs s1:τ, and we neglect possible uncertainty arising from s1:τ, thereby allowing us to ignore the corresponding density p(s1:τ) in Eq. (7). The likelihood p(y1:τ|x1:τ,ϕ,s1:τ) determines how likely the observational data are, given a sequence of states produced by the forward model, while the distribution p(x1:τ|ϕ,s1:τ), sometimes referred to as the “process model”, relates a particular set of parameters ϕ to a sequence of states x1:τ. In our approach, we focus on probabilistically modeling parametric uncertainty induced by the prior p(ϕ)=i=1mpi(ϕi) over m model parameters which encode a priori site-specific knowledge (e.g., soil characteristics and surface conditions) where available. We neglect uncertainty in the process model, thereby assuming the model states to be a deterministic function of the parameters. Additional theoretical discussion of the inverse modeling method can be found in Appendix B3.1.

3.4.2 Approximate inference via ensemble Kalman sampling

Evaluation of the likelihood p(y1:τ|x1:τ,ϕ,s1:τ) requires a forward simulation of the heat transfer model described in Sect. 3.3 given the parameters ϕ which incur a substantial computational cost. To mitigate this issue, we employ a recently developed, probabilistic variant of the widely used ensemble Kalman inversion (EKI) (Iglesias et al.2013) algorithm called ensemble Kalman sampling (EKS) (Garbuno-Inigo et al.2020) to efficiently produce approximate samples from the posterior over model parameters, p(ϕ|y1:τ,s1:τ). EKS assumes that the observed borehole temperatures can be represented as

(8) T obs = ( g T f ) ( ϕ ) + η ,

where gT is the forward map from the model states produced by f to comparable temperature observations, and ηN(0,ΣT) is observation noise with zero mean and assumed covariance ΣT. Additional inputs to f, such as forcing data, are assumed to be independent of ϕ and are here suppressed for brevity. EKS requires the prior distribution over parameters to be a multivariate Gaussian, so we must first construct a bijective function γ:ΦΨ that maps the m-dimensional parameters ϕ∈Φ to values ψ∈Ψ with unbounded support on the real line. The approximate prior for EKS is then specified as ψp(ψ)=N(μψ,Σψ). The key insight behind the EKS algorithm is that an ensemble of parameter values {ψ(i)} sampled from an initial density (i.e., the prior) can be transformed into samples from the (approximate) posterior distribution by treating them as a stochastic system of interacting particles. Their dynamics are then governed by the overdamped Langevin equation:

(9) ψ ( i ) t = - L ( ψ ( i ) ) + 2 R ( t ) ,

where ℒ(ψ) is the log-likelihood of the unconstrained parameters ψ with respect to the observations and R(t) is an m-dimensional Brownian motion. It can be shown that this interacting particle system then converges to the posterior density over an infinite time horizon, where “time” here refers to that of the particle diffusion rather than physical time in the forward model. We refer the reader to Garbuno-Inigo et al. (2020) for further details about the EKS algorithm and iteration procedure.

The main advantage of EKS over other commonly used numerical sampling methods, such as Markov Chain Monte Carlo (MCMC), is that it is a parallelizable “ensemble” or particle-based method which requires a relatively small number of iterations to converge (typically fewer than 30). In contrast, most MCMC methods require thousands of sequential (and therefore non-parallelizable) iterations. This makes EKS more attractive for use cases such as the inverse modeling problem considered in this work where the computational cost of likelihood is relatively high. EKS also avoids the need to design complex, multi-stage model reduction and/or emulation procedures such as those used for Bayesian evidential learning (Hermans et al.2018). This comes at the cost, however, of finding only a Gaussian approximation to the “true” posterior.

We apply the implementation of EKS published by the authors (Constantinou et al.2022) to our heat conduction model using borehole data obtained from each of the four sites described in Sect. 2. We use an ensemble of size N=512, and we initialize the ensemble by sampling N parameter vectors ϕ from the prior p(ϕ). We iterate EKS for at most 30 iterations or until two units of particle diffusion “time” are reached. Further details on the problem setup for EKS are given in Appendix B3.2.

Figure 4Trend analysis of observed mean annual temperature trends in air vs. soil at all sites. Permafrost temperatures are from borehole measurements as close as possible to the depths of zero annual amplitude: 20.75 m (Samoylov), 20.0 m (Barrow), 9.0 m (Bayelva), and 9.82 m (Parson's Lake). The dots represent observed mean annual temperatures at each location. For air temperature, the adjusted ERA5-Land temperature reanalysis record is used where observations are not available. Trend lines and their corresponding 95 % central credible intervals are obtained via numerical sampling from a robust trend model.


4 Results

4.1 Observed permafrost and air temperature trends

The results of the trend analysis over the 20 year time period (2000–2020) described in Sect. 3.2 are shown in Fig. 4. The left panel shows the observed air and permafrost temperatures over time alongside the trend lines sampled from the posterior. The solid line represents the posterior mean, while the shaded region shows the 95 % central credible interval (CCI), which represents the interval between the 2.5 % and 97.5 % percentiles of the posterior samples. The boxplots in the right panel show summaries of the posterior distribution over slopes of the trend line, which can be more naturally interpreted as the average annual change in temperature over the observed time period. The whiskers of the boxplot span the range [q1-1.5R,q3+1.5R], where q1 and q3 are the first and third quartiles, respectively, and R is the interquartile range. All probabilities given in this section refer to posterior probabilities calculated from MCMC samples, with Monte Carlo adjusted standard error (MCSE) strictly less than 2×10-3 in all cases.

Our analysis shows very strong evidence of increasing air temperatures (i.e., Pr(μ1>0)>0.99 for all sites), with the slope of the trend line μ1 ranging from 0.09 K yr−1 (95 % CCI 0.00 to 0.18 K yr−1) at Parson's Lake to 0.11 K yr−1 (0.04 to 0.18 K yr−1) at Samoylov Island. In comparison, mean annual permafrost temperatures near the depth of zero annual amplitude show similarly strong evidence of warming at three of the four sites: temperatures rose by approximately 0.16 K yr−1 (0.14 to 0.18 K yr−1) on Samoylov Island, 0.06 K yr−1 (0.05 to 0.06 K yr−1) at Barrow, and 0.06 K yr−1 at Parson's Lake (0.02 to 0.11 K yr−1). For Bayelva, the observed deep permafrost temperatures show a far milder warming trend of approximately 0.02 K yr−1 (0.02 to 0.05 K yr−1) and are considerably less conclusive (Pr(μ1>0)0.85).

Next, we investigate the observed difference between soil and air temperature trends across sites. We find that soil temperature trends are generally lower than air temperature trends at both Barrow (Pr(μ1(soil)<μ1(air))0.90) and Bayelva (Pr(μ1(soil)<μ1(air))0.98). This relationship is less evident at the Parson's Lake borehole site (Pr(μ1(soil)<μ1(air))0.68), with ground temperatures warming at a comparable rate to air temperatures. In stark contrast to the other sites, borehole measurements at Samoylov show deep soil temperatures warming on average 0.05 K yr−1 (35 %) faster than air temperatures with fairly high confidence (Pr(μ1(soil)>μ1(air))0.92).

Figure 5(a) Modeled vs. observed mean annual ground temperatures at sensors which are as close as possible to the depth of zero annual amplitude: 20.75 m (Samoylov), 20.0 m (Barrow), 9.0 m (Bayelva), and 9.82 m (Parson's Lake). The dots indicate observed annual borehole temperatures, and the lines represent the mean of the ensemble predictions. The shaded regions cover 95 % of the ensemble predictions. (b–e) Minimum and maximum annual temperature profiles at 2 m and below for the years 2020 (Samoylov, Bayelva, Parson's Lake) and 2019 (Barrow), observed (dots) vs. model ensemble predictions, with the white dotted line indicating the median over the ensemble.


4.2 Ensemble inversion

Modeled vs. observed annual ground temperatures over the 20 year simulation period are shown in Fig. 5. The predicted temperatures from each ensemble member represent an approximate sample from the posterior predictive distribution. The solid line in Fig. 5a is the ensemble mean temperature prediction, while the shaded regions show the 95 % CCI over all ensemble members, i.e., the interval into which 95 % of the predicted values for each ensemble member fall. Figure 5b–e show the observed vs. ensemble minimum and maximum annual temperatures for the last year in the simulation period where observations are available. We can see in Fig. 5 that, for both Samoylov and Parson's Lake, the ensemble predicted temperatures fit the observed borehole temperatures relatively well, although the 2020 annual min–max ground temperature curve for Samoylov (Fig. 5b) indicates a slight bias in the min–max temperature range for the upper part of the profile affected by seasonal variation. The ensemble temperatures for Bayelva also appear reasonable, though there is a clear warm bias in the winter apparent in the temperature profile. The corresponding min–max temperature curves (Fig. 5d) reveal a substantial amount of ensemble variation in the annual minimum and maximum temperature despite fairly modest variation in the ensemble mean at 9 m. This is also apparent in the 2020 ensemble temperature profile for Parson's Lake (Fig. 5e). There is much more variation in the deeper part of the soil profile (i.e., below 10 m) for both of these sites in comparison to Samoylov and Barrow (Fig. 5b–e), likely due to the lack of deeper borehole measurements available to constrain the initial temperature profile.

Figure 6Joint densities of modeled mean annual change in latent heat (a) and active layer thickness (b) vs. mean annual change in ground temperature for all sites. Each point represents a single member of the fitted ensemble. Energy content and temperature are integrated over all grid cells in the upper 10 m of the soil profile. The density plots on the top and right panels of each plot show the marginal densities. In both plots, a small number of points from Bayelva exceed the upper y-axis limit and thus are not shown.


Figure 6 shows the joint and marginal distributions of average annual change in modeled latent heat (Fig. 6a) and ALT (Fig. 6b) with average integrated ground temperature, i.e, the weighted average of grid cell temperatures, for each of the four sites. Here, average annual change is defined as the slope of the least squares trend line fitted over the entire simulation period (2000–2020) for each ensemble member. We integrate the predicted temperatures and energy contents over the upper 10 m of the soil profile. The choice of 10 m is arbitrary but is motivated by the availability of measurements from Bayelva and Parson's Lake. Note that, since the model assumes constant hydrological conditions, changes in latent heat can be directly attributed to phase change. The density plots on the top and right panels of Fig. 6 show the ensemble marginal densities of temperature and latent heat change, respectively. Each point in the plot represents the slope of the least squares trend line for one ensemble member after fitting with EKS. The lines in the plot show the linear relationship between changes in temperature and latent heat across all ensemble members. Table 2 shows the correlation and regression coefficients associated with each linear fit. Stronger positive correlations between temperature and latent heat trends imply higher sensitivity of permafrost thaw to changes in sensible heat uptake.

Table 2Summary statistics (Pearson correlation, slope of the least squares line, and standard deviations) for the linear relationships shown in Fig. 6 with respect to both latent heat vs. temperature (a) and ALT vs. temperature (b). Means and standard errors are bootstrap estimates using 10 000 resamples from the ensemble outputs.

Download Print Version | Download XLSX

We can see in Fig. 6a that all three sites with silty or clay-like soils (i.e., Barrow, Bayelva, and Parson's Lake) exhibit a much stronger positive correlation between change in latent heat and temperature. Samoylov, in contrast, shows little to no correlation despite the comparatively rapid increase in ground temperature. Figure 6b shows a clear difference in how changes in ALT correlate with changes in temperature between the two warmer sites, Parson's Lake and Bayelva, versus the two colder sites, Samoylov and Barrow. The simulations for both Parson's Lake and Bayelva show a stronger positive correlation between increasing ALT and increasing temperature, whereas for the colder sites this correlation is weaker (Barrow) or virtually non-existent (Samoylov). Both warmer sites also show more overall variability across the ensemble in both latent heat and ALT trends, with standard deviations of the trend slopes more than double those of the cold sites.

Figure 7Total change in energy partitions (upper 10 m of soil profile) since the beginning of the simulation period. Energy is partitioned into three categories: latent heat in frozen grid cells (i.e., cells with maximum annual temperature Tmax< 0 C), latent heat in the active layer (Tmax 0 C), and sensible heat. Solid lines show the median energy change, while the shaded regions show the 95 % CCI over the ensemble. Changes in heat capacity due to phase transitions in pore water and ice are excluded here from the sensible heat component in order to better isolate the relative contributions of temperature change vs. phase change.


Lastly, Fig. 7 shows the long-term change in the mean annual energy contents of the upper 10 m of the soil profile over the full simulation period (2000–2020) for each site. The energy storage is partitioned into three parts: latent heat stored in permafrost, i.e., grid cells that have a maximum annual temperature below 0 C; latent heat stored in the active layer; and sensible heat due to temperature change. The solid lines show the median energy change over all fitted ensemble members, and the shaded regions indicate the corresponding 95 % CCI. The categorization of a grid cell as being in the permafrost or active layer for the latent heat partitions is not constant; i.e., as the active layer grows or shrinks over time, the categorization of each grid cell based on maximum annual temperature changes as well. This means that the change in latent heat for each partition is attributable not only to the thawing and freezing of water but also to the migration of grid cells between the two partitions.

Note that the sensible heat as defined in Eq. (1) depends on the liquid water content via the heat capacity (see Eq. B6 in the Appendix). For the purpose of more clearly visualizing the relative contributions of temperature (sensible) vs. latent heat change in Fig. 7, we employ a non-conservative formulation of the sensible heat: Hs(T)=T×ci, where ci is the volumetric heat capacity of ice, thereby effectively removing the dependence of sensible heat on the unfrozen water content. This alternative formulation allows us to more clearly see the relative contributions of temperature change and latent heat to the long-term evolution of the thermal state.

5 Discussion

5.1 Estimating trends in observed permafrost and air temperatures

Numerous recent studies have highlighted warming trends observed in borehole sites all across the Arctic (Smith et al.2022; Chen et al.2021; Biskaborn et al.2019; Romanovsky et al.2007). The results of our trend analysis are broadly consistent with Biskaborn et al. (2019), who reported average decadal changes in permafrost temperatures of 0.39 ± 0.15 K dec−1 for 53 GTN-P boreholes across the continuous permafrost zone between 2007 and 2016, which translates into an average annual change of roughly 0.02 to 0.05 K yr−1. Observed warming trends at the Samoylov Island borehole site in the Lena River delta (Siberia) over the longer period of 2007–2020 are unusually extreme (also noted by Biskaborn et al. (2019) for the same borehole), with a mean annual change of 0.16 K yr−1, while observed changes at the Bayelva site near Ny-Ålesund (Svalbard) fall towards the lower end of the expected range with a mean annual change of 0.02 K yr−1. For air temperature, Biskaborn et al. (2019) report average changes of 0.86 ± 0.84 K dec−1 (0.02 to 0.17 K yr−1) across the continuous permafrost zone from 2007–2016 using spatially interpolated 2 m air temperatures from reanalysis products. This is consistent with the results of our analysis which show mean air temperature trends ranging from 0.09 K yr−1 at Parson's Lake to 0.11 K yr−1 on Samoylov Island. This discrepancy in observed permafrost temperature trends between these two sites, despite similar changes in air temperature, indicates considerable uncertainty in how permafrost is responding to the changing climate that is likely attributable to other factors. For example, thicker and/or less dense snow cover can accelerate permafrost warming by insulating the ground against rapid drops in air temperature characteristic to autumn and early winter, thereby delaying the refreezing of the active layer (Park et al.2015). Additionally, soil thermal properties such as the bulk conductivity and freezing characteristics due to soil texture can also play a significant role in modulating the effects of surface temperature changes (Nicolsky and Romanovsky2018; Romanovsky and Osterkamp2000; Riseborough1990). Both of these factors, among others, can significantly affect energy uptake in the subsurface and, ultimately, the current and future thermal state of permafrost in Arctic regions (Smith et al.2022; Langer et al.2022).

The results of this trend analysis motivate our inverse modeling study by demonstrating clear, localized differences and substantial uncertainty in how the permafrost thermal regime at these four sites is responding to long-term changes in air temperature. We argue that this can be at least partially attributed to factors affecting the uptake of latent heat in the subsurface such as soil freezing characteristics as well as historical climatology. We discuss in the following sections how our inverse modeling results suggest that both of these factors are major sources of uncertainty in making inferences about the changing thermal state of permafrost.

5.2 Variability and biases in modeled permafrost temperatures

There is substantial variability in modeled temperatures over the simulation period (2000–2020) even after using EKS to calibrate the model ensemble to borehole measurements (Fig. 5). We see, for example, significant differences in the overall range of predicted temperatures between sites. This can be partially explained by both differences in air temperature variability (Fig. 4) as well as in the depths at which temperature observations are available; e.g., the deepest sensors at the Bayelva and Parson's Lake sites are at 9 m and 9.82 m, respectively, whereas both Barrow and Samoylov have much deeper measurements available near the depth of zero annual amplitude where there is little to no impact from seasonal variation. We present additional experiments exploring the sensitivity of our inversion method to available sensor depths in the Supplement (Sect. S1 in the Supplement). We find that the inclusion or exclusion of measurements at various depths has a fairly predictable effect on the results of the analysis; i.e., the inclusion of data at additional sensor depths generally improves the fidelity of the predicted temperatures at these depths and reduces the associated uncertainty. However, the results suggest that constraining the model using only deep measurements, as we do for the Barrow site, may lead to overestimation of warming at depths closer to the surface.

The general agreement between the modeled and observed mean annual temperatures in Fig. 5 suggests that the fitted model ensemble can nonetheless provide a plausible estimate of the permafrost's thermal state at all four sites, subject to some uncertainties. There are, however, some clear discrepancies between modeled and observed permafrost temperatures which should be considered, particularly in the final year temperature profiles for Samoylov and Bayelva. For Samoylov, the ensemble annual temperature range of the permafrost in the upper 10 m of the ground is slightly too narrow, showing a warm bias in the winter (0.6 ± 0.3 K averaged over the profile) and a cold bias (-0.7±0.2K) in the summer for the last simulation year of 2020 (Fig. 5b). This may be due to the thermal conductivities of the upper layers being too low, thereby over-insulating the lower layers. Furthermore, while many of the soil properties are varied in the ensemble, the stratigraphy itself (i.e., the number and thickness of layers) is not, which could also be responsible for such biases. The wintertime warm bias for Bayelva (1.5 ± 0.5 K) may also be due to our assumption of static hydrology, which fails to capture the effects of seasonal changes in soil moisture on the subsurface thermal regime. The release of latent heat during freezing slows the propagation of cold surface temperatures, thereby delaying the refreezing of the active layer in early winter (Romanovsky and Osterkamp, 2000). When unfrozen water is removed from the active layer due to drainage or evapotranspiration, the latent heat stored in this water is removed as well. The result is that cold temperatures can propagate faster in the wintertime since a higher fraction of this energy is diffused as sensible heat. This would therefore serve as a possible explanation for the warm bias in the model at this site. We can see quite clearly from the volumetric water content measurements shown in Fig. 3 of Boike et al. (2018a) that soil saturation levels in the active layer vary substantially (10 % to 20 %) both within and between years. It has also been well established by previous studies that hydrological conditions have a substantial effect on the subsurface thermal regime on Svalbard (Westermann et al.2011; Roth and Boike2001). Accounting for water fluxes may therefore be necessary to more accurately reproduce thermal conditions at Bayelva and other sites with similar characteristics. Nevertheless, we believe the results of our inversion procedure agree well enough with temperature observations to provide insight into how the thermal state of permafrost at each site may have evolved over the last 2 decades.

5.3 Linking temperature trends to changes in permafrost thaw

The clear differences between the joint distributions for each site in Fig. 6 have substantial implications for the interpretation of permafrost temperature trends. Firstly, it is evident that the freezing characteristics of the soil, represented in our inversion method through variation of the soil freezing characteristic curve parameters, have a substantial impact on the relationship between trends in temperature and latent heat. This is made most apparent in Fig. 6a by the presence of a strong linear relationship between latent heat change and temperature change at Barrow (18.9 ± 0.8 MJ K−1, Table 2) but not at Samoylov (2.1 ± 0.9 MJ K−1). Both sites are in many ways similar; they both have historically very cold permafrost (below 8 C near the depth of annual amplitude), an insulating organic layer at the surface, and similar amounts of annual snow cover (Table 1). One of their key differences, however, is that the soils in the Barrow region tend to have higher silt content (Romanovsky and Osterkamp2000) which exhibits freezing characteristics that lead to larger amounts of unfrozen water at subzero temperatures. The result is that, as permafrost close to the active layer warms, a larger amount of energy is consumed as latent heat, thereby slowing the rate of temperature change near the freezing front and effectively insulating the colder permafrost below (Nicolsky and Romanovsky2018). This is also seen clearly in Fig. 7, where the estimated change in latent heat storage of permafrost for the Barrow simulations clearly dwarf that of Samoylov by up to a factor of 10.

The joint density of active layer thickness trends with temperature trends in Fig. 6b tells a similar story. As with latent heat, simulations for the three sites with silty or clay-like soils in the permafrost layers show the strongest positive correlation with ALT, while those for Samoylov show little to none. However, the correlation and effect size for the Barrow simulations (0.26 ± 0.04 m K−1 and 0.03 ± 0.01 m K−1, respectively, Table 2) are notably lower than those for the two warmer sites, Bayelva (0.35 ± 0.03 m K−1 and 0.13 ± 0.02 m K−1) and Parson's Lake (0.46 ± 0.06 m K−1 and 0.06 ± 0.01 m K−1). This is likely due to the colder, deep permafrost at Barrow acting as a heat sink and thereby offsetting the effects of warming closer to the surface.

The vast majority of ensemble members for both Samoylov and Barrow show clear warming trends, whereas the marginal distributions of temperature change for Bayelva and Parson's Lake are more diffuse and tend to be closer to zero. This is broadly consistent with the hypothesis of Riseborough (1990) as well as the more recent findings of Nitzbon et al. (2023), i.e., that regions with warmer permafrost may show less detectable subsurface warming while actually being at higher risk for degradation. The two warmer sites (particularly Bayelva) also tend to show similar or larger changes in latent heat content and ALT even with relatively little change in temperature; this is likely due to a larger portion of energy being consumed as latent instead of sensible heat. This hypothesis is further supported by Fig. 7, which shows sensible heat as the clear dominant component of long-term energy change at Samoylov and, to a lesser extent, also at Barrow. We discuss this further in Sect. 5.4.

The wider range of modeled ALT and latent heat change at both Bayelva and Parson's Lake (Figs. 6 and 7) indicates higher uncertainty in modeling the thermal dynamics of warmer permafrost. This can be attributed to higher sensitivity to soil properties, initial conditions, and changes in surface conditions and soil water content. This sensitivity is in large part due to the nonlinear effects of the freeze curve on the thermal dynamics, especially in deep permafrost (i.e., near the depth of zero annual amplitude) where thermal gradients are smaller and heat diffuses more slowly.

A small number of ensemble members for both of these warmer sites lie in the upper left quadrants of Fig. 6a and b, i.e., where there is a positive trend in latent heat or ALT despite average ground temperatures actually declining over the simulated time period. Similarly, some ensemble members for these two sites fall into the lower two quadrants, indicating a decrease in latent heat or ALT with some parameter settings. This is in contrast to the two colder sites, Samoylov and Barrow, where nearly all ensemble members show positive trends in temperature, latent heat, and ALT. This underscores two features of warmer permafrost: (i) that observed temperature trends (or lack thereof) should be interpreted with caution, as there is substantial uncertainty inherent in associated changes in latent heat, and (ii) that the sensitivity of warmer permafrost to climatic change is highly dependent on soil thermal properties, most notably the soil freezing characteristics. Soils which retain more unfrozen water at temperatures below 0 C may delay thawing and refreezing due to the nonlinear effects of latent heat (Nicolsky and Romanovsky2018).

5.4 Latent heat as a subsurface energy sink

While it is clear that Arctic permafrost is warming (Smith et al.2022; Chen et al.2021; Biskaborn et al.2019; Smith et al.2010), the latent heat effect acts on a larger scale as both an energy source and sink (Nitzbon et al.2023; Cuesta-Valero et al.2023), which can have the effect of masking long-term temperature changes in the climate signal (Nicolsky and Romanovsky2018; Riseborough1990), thereby introducing considerable uncertainty into the long-term implications of warming across the Arctic. Figure 7 visualizes this uncertainty by comparing long-term changes in the partitioned subsurface energy budget across sites. There is a clear contrast between the two colder sites, Samoylov and Barrow, in the upper panels versus the two warmer sites, Bayelva and Parson's Lake, in the lower panels. For the colder sites, sensible heat change plays a more dominant role with both the median and full 95 % CCI of the ensemble showing a positive trend for most of the simulation period. The decrease in sensible heat at Barrow in the last simulation year, 2020, is likely due to the unusually large drop in air temperature which can be seen in Fig. 4. In contrast, both warmer sites show latent heat change playing a larger role.

Long-term change in latent heat uptake also differs between the colder and warmer sites. For both Samoylov and Barrow, a vast majority of ensemble members show an increase in latent heat both in the permafrost partition (i.e., grid cells where maximum annual temperatures are below zero) as well as in the active layer. This points to both a deepening active layer (consistent with Fig. 6) as well as to warming permafrost. For sites like Barrow with silty soil, the notably larger increase in permafrost latent heat storage also implies a transition of the permafrost into a kind of “mushy” (partially frozen) state as temperatures increase into the critical region of the soil freezing characteristic curve, where liquid water is present at temperatures below the nominal melting point.

Both Bayelva and Parson's Lake show a much larger range of possible trajectories for latent heat storage over the ensemble. Notably, however, the 95 % credible range for both sites (especially Bayelva) spans a wide range of possible trajectories for latent heat storage in the permafrost. Latent heat in the active layer, however, appears to be likely increasing at both sites. This may be attributable to expansion of the active layer in combination with delayed sensible warming of permafrost due to the latent heat effect. Crucially, both warmer sites also consist of a mixture of more silty and clay-like soils which are characterized by larger amounts of unfrozen water content at temperatures well below the melting point (Fig. 2a). One interesting difference even between the two warmer sites is the relative contribution of latent heat change in the permafrost layer (light blue in Fig. 7) versus sensible heat (light red). The larger variation in this latent heat component at Bayelva appears to be due to two factors: (i) warmer permafrost temperatures in the upper 10 m (Fig. 5) and (ii) differences in the soil stratigraphies; e.g., our model configuration for Parson's Lake has on average 10 % lower porosity from 1.75 to 10 m depth which translates into less overall latent heat by volume. Furthermore, the more heterogeneous mixture of soil textures (and therefore freezing characteristics) in the Parson's Lake stratigraphy, as well as the higher variability in air temperature (Fig. 4), results in a higher overall contribution of sensible heat to the subsurface energy balance. Our analysis here highlights how the nonlinear effect of the soil freezing characteristics on permafrost thermal dynamics limits the reliability of temperature measurements in understanding the impact of climate warming on the thermal state of permafrost.

5.5 Limitations

While we are confident that the ensemble inversion method and results presented in this study provide valuable insight into the relationship between air temperature trends, subsurface temperature trends, and permafrost thaw, there are several limitations of our methodology which are important to address.

  1. Neglected subsurface processes. The transient heat conduction model described in Sect. 3.3 is capable of accurately simulating heat conduction with phase change at a relatively high spatial resolution. However, it neglects several processes such as water infiltration and percolation, excess ice melt and subsidence, and lateral exchange of heat and water in the subsurface. In particular, the assumption that soil water content remains constant over the simulation period precludes factoring in the effects of long-term wetting or drying on the thermal properties (and therefore the thermal regime) of the ground closer to the surface. This could potentially alter the thermal dynamics of the ground in two ways: firstly, unfrozen water which is transported out of the pore space and not replaced can be seen as effectively removing energy from the system, thereby altering the energy balance of the ground. Secondly, drier soils with higher air content will necessarily be less conductive, thereby becoming an insulator for deeper layers against further warming and cooling. However, the assumption of static hydrological conditions allows us to more easily compare changes in the partitioning of latent vs. sensible heat between sites (Fig. 7). We further note that for Samoylov, which is typically waterlogged, no substantial wetting or drying has been reported during the simulation period (Boike et al.2019).

  2. Simplified parameterization of surface processes. Our model does not explicitly represent surface processes such as radiative and turbulent heat fluxes, water bodies, vegetation, and snow dynamics but represents them in bulk via n factors (see Sect. 3.3.2). Thus, inter-annual variability in surface conditions is not effectively represented, and meteorological forcings other than air temperature are not taken into account, thereby neglecting potentially valuable sources of information that could serve as additional constraints on the inverse problem. We might expect, for example, some of the inter-annual fluctuations to be explained by changes in snow or cloud cover. The inclusion of additional surface processes, however, also brings with it a whole range of additional uncertainties from both the forcing data as well as from the additional parameterizations.

  3. Selection and calibration of prior distributions. The model parameter priors employed in this work are loosely derived from both auxiliary data sources as well as field work and published values in the literature (see Appendix B3.3 for further discussion). However, the general lack of precise error and uncertainty estimates for some model parameters (e.g., soil properties) makes the construction of prior distributions difficult. It is evident from the marginal posterior plots (Figs. A2A5) that many parameters show minimal change from their priors. This may be the result of overly narrow prior densities or may simply indicate that annual temperature measurements do not provide sufficient information to update them. Furthermore, it is very likely that many parameters are, in reality, highly correlated in the joint distribution, which we are currently unable to account for in the prior. Prior distributions which are more informative or more realistic will help the inversion method converge faster and improve the reliability of the resulting posterior samples in making inferences about the system under study. We nevertheless emphasize that even imperfect (but plausible) prior distributions still do a better job of accounting for uncertainty than the point estimates of parameter values typically provided by non-probabilistic calibration methods.

  4. Statistical limitations. While EKS provides a computationally efficient alternative to MCMC for drawing samples from the posterior, it is nevertheless an approximation which provides a theoretical guarantee of convergence to the posterior measure only over an infinite “time” horizon (i.e., iterations of the algorithm) (Garbuno-Inigo et al.2020). Furthermore, the empirical results presented by the authors of EKS indicate that the method underestimates posterior variance (and thus uncertainty) on a finite time horizon unless the ensemble is very large. More recent work by Cleary et al. (2021) has attempted to circumvent this problem by performing exact posterior inference on emulated predictions produced by a surrogate model. Additionally, EKS requires a priori specification of the sampling distribution covariance ΣT, which is not, in general, known. A fully Bayesian treatment of the inversion problem would include these ΣT (or equivalent) parameters in the posterior distribution. This would have the benefit of producing posterior samples where the predictive distribution is well calibrated; i.e., the 95 % prediction interval should actually cover approximately 95 % of the observations in the calibration period, therefore providing a built-in measure of the degree to which the model is capable of explaining the observations.

These limitations provide a wide range of directions for future work; e.g., points 1 and 2 could be addressed simply by incorporating additional physical processes into the model, such as those described in Westermann et al. (2023). However, this typically comes at a heavy computational cost, limiting the feasibility of running large ensembles over longer time periods spanning multiple decades, centuries, or millennia. We are optimistic that recent advancements in simulation-based inference (Cranmer et al.2020) and/or machine-learning-based approximations of the likelihood (Thibaut et al.2022; Hermans et al.2018) may be able to address this challenge, as well as those described by issues 4 and 5. We also see potential in applying an inversion method similar to the one outlined in this work at a larger, circum-Arctic scale, perhaps by leveraging more efficient formulations of the heat conduction model such as those proposed by Tubini et al. (2021) or Langer et al. (2022). The greatest challenge here is the limited availability of high quality borehole data at a global scale, despite the ongoing efforts of collaborative data collection and publication projects such as GTN-P (Biskaborn et al.2015).

6 Conclusions

In this work, we analyzed trends in annual air and permafrost temperatures at four different sites in the Arctic continuous permafrost zone. We found that long-term trends in permafrost temperatures over the study period (2000–2020) varied between the four sites (0.02 to 0.16 K yr−1) more than trends in air temperatures (0.09 to 0.11 K yr−1). In order to investigate what permafrost temperatures can tell us about the thermal state of permafrost, we then applied a recently developed method for ensemble inverse modeling, ensemble Kalman sampling (Garbuno-Inigo et al.2020), to constrain a transient heat conduction model using long-term borehole temperature measurements from each site. We demonstrated that the response of the subsurface thermal regime to warming at the surface varies substantially between sites based on their historical climatology (i.e., colder versus warmer conditions) as well as the thermal properties of the soil. We found that the degree to which trends in permafrost temperatures correlate with changes in latent heat or active layer thickness depends heavily on both historical climatology and the soil freezing characteristics. In particular, the sensitivity of modeled active layer thickness to changes in temperature were higher for the two sites with warmer permafrost, 6 ± 1 cm K−1 at Parson's Lake in northwestern Canada and 13 ± 2 cm K−1 at Bayelva, Svalbard, versus 3 ± 1 cm K−1 at Utqiaġvik (Barrow) in northern Alaska and a negligible 1 ± 1 cm K−1 at Samoylov Island in northeastern Siberia. Furthermore, we showed that there is substantial uncertainty inherent in relating ground temperature trends to permafrost thaw in areas with warmer permafrost or silty or clay-like soils where the nonlinear effects of latent heat play a much larger role. Our results suggest that currently observed permafrost temperature trends at these sites can be plausibly mapped to a wide range of different thermal states. This is largely due to the lack of information in the temperature signal of permafrost within the critical temperature range where phase change of water occurs. Our findings demonstrate the importance of considering regional differences in soil characteristics and climatology when comparing permafrost temperature trends across sites. Our study also highlights the importance of deep borehole measurements (i.e., deeper than 10 m) in constraining hydrothermal models of permafrost and providing insights into the past and present thermal state of permafrost in the Arctic. Future work will focus on the broader application of probabilistic inversion methods, like the one used in this work, both to more sophisticated permafrost models and at larger spatiotemporal scales.

Appendix A: Additional results

A1 Posterior distributions of wintertime n factors indicate changes in snow cover

The posterior distributions of the wintertime n factors are shown in Fig. A1. Each group of boxplots represents the fitted value of the n factor knot at that year. As described in Sect. 3.3.2, the actual n factor values are interpolated linearly between knots, forming a piecewise linear function over the simulation period. Values for the year 2000 are also applied for all years in the spin-up period (1979–1999). The inter-annual differences between the knots indicate that, in addition to the long-term trend in air temperature, changes in surface conditions are necessary in order to accurately reproduce temperature trends in the borehole data.

Decreasing wintertime n factors for Samoylov between 2005–2020 can be physically interpreted as an increase in surface insulation, typically explained by an increase in snowfall or changes in the timing of the snow covered season. The extent to which long-term snow depth measurements published in Boike et al. (2019) support this hypothesis is unclear, and further validating it will likely require further investigation with models that take into account snow dynamics over shorter timescales. We note, however, that the long-term pattern produced by the n factors serves as a possible explanation for the rapid rate at which deep ground temperatures appear to be warming at the Samoylov Island site (Fig. 4). As discussed in Boike et al. (2019), a new research station was constructed within the immediate vicinity of the borehole site in 2013. The station includes a water delivery system which may act as a wind shield, thereby affecting seasonal snow cover patterns of the area around the borehole. Visual inspection in April of 2016 indicated that snow cover did indeed appear higher since construction of the water supply system (Boike et al.2019). Thus, there is reason to believe that snow cover has increased at this site, though the cause of this increase remains unclear.

Figure A1Boxplots showing posterior distributions of wintertime n factors (i.e., applied when Tair(t)≤0) after fitting with EKS.


For the site at Bayelva, Svalbard, the fitted wintertime n factors show the opposite effect, i.e., increasing from roughly 0.4 in 2000 to roughly 0.6 in 2015. This can be interpreted as a reduction in surface insulation, suggesting decreasing snow cover or possibly a longer snow-free season. Long-term snow observations presented by Boike et al. (2018a) verify this hypothesis. The fact that fitted n factor distributions recover this pattern suggests that changes in snow cover are necessary in order to reconcile observed borehole temperatures with warming air temperatures. The higher posterior variance of the Bayelva n factors vs. Samoylov also indicates more uncertainty in the strength of the effect likely due to larger relative contribution of latent heat limiting the information gain from temperature measurements. It is also worth noting that the n factor distributions for Barrow and Parson's Lake show less of a clear temporal pattern, indicating that the observed borehole temperatures can plausibly be explained by air temperature warming alone.

Table A1Names and descriptions of parameters considered in the ensemble inversion procedure. Some parameters appear multiple times across years or layers as indicated.

Download Print Version | Download XLSX

A2 Parameter prior vs. posterior distributions

Histograms of posterior samples (i.e., ensemble members) for all parameters and for each site are presented in Figs. A2 to A5. Parameter names are summarized in Table A1 and generally follow the pattern[.sub], where name refers to the model parameter being determined and sub refers to the transformed “sub-parameter” where applicable. We use a different naming convention for the forcing and initial condition parameters in order to clearly distinguish them. Note that the set of varied parameters differs slightly between sites due to differences in the stratigraphies; e.g., the excess ground ice parameter is not considered in sites which are assumed (correctly or not) to have no excess ice within the soil column. We can see that, in most cases, the prior and the posterior generally agree, with the only exception being the 2010 and 2015 wintertime n factors for Samoylov. There is clearly a larger reduction in uncertainty of the n factor parameters for both Samoylov and Barrow which is consistent with our findings that sensible heat change dominates at both sites and thus provides a clear signal for changes in the air temperature forcing. Both Parson's Lake and Bayelva generally show minimal change to the prior distribution, likely due to the dominant contribution of latent heat which limits the information gain from temperature measurements.

Figure A2Fitted ensemble posterior parameter distributions – Samoylov Island.


Figure A3Fitted ensemble posterior parameter distributions – Barrow.


Figure A4Fitted ensemble posterior parameter distributions – Bayelva.


Figure A5Fitted ensemble posterior parameter distributions – Parson's Lake.


Appendix B: Method details

B1 Trend analysis

The annual temperature trend model is specified as


where T(μ,σ,ν) is the scaled and shifted Student's t distribution with ν degrees of freedom, and t is the zero-centered year. Samples from the posterior distribution p(μ1,μ0,σ,ν|Tobs,t) are obtained via the widely used “No-U-Turn” numerical sampler (NUTS) (Hoffman and Gelman2014). We combine 500 independent samples from eight different chains with 2500 discarded adaptation samples to get a total of n=4000 samples for each trend model.

B2 Transient heat conduction model

B2.1 Soil parameterization

The soil constituent fractions by volume are parameterized using the four characteristic fractions: excess ice (ξxic), natural porosity (ξpor), saturation (ξsat), and organic solid fraction (ξorg). The volumetric soil constituent fractions for excess ice (θx), pore water and ice (θwi=θw+θi), air (θa), and mineral and organic content θm,θo can be computed according to the following linear system:


where θx+θwi+θa+θo+θm=1. Mean stratigraphy configurations (i.e., expected value of the prior) used for each site are given in Tables B1 to B4. Note that, in some cases, adjacent layers have identical mean values but are varied independently in the ensemble. Both Bayelva and Parson's Lake also have different freeze curve priors for some layers based on differences in observed soil texture; this is discussed further in Appendix B3.3.

Table B1Stratigraphy configuration for Samoylov Island borehole simulations.

Download Print Version | Download XLSX

Table B2Stratigraphy configuration for Bayelva borehole simulations.

Download Print Version | Download XLSX

Table B3Stratigraphy configuration for Barrow borehole simulations.

Download Print Version | Download XLSX

Table B4Stratigraphy configuration for Parson's Lake borehole simulations.

Download Print Version | Download XLSX

We use the soil freezing characteristic of Dall'Amico et al. (2011), which is based on the soil water retention equation of van Genuchten (1980). The temperature dependent liquid water matric potential is expressed as

(B1) ψ ( T ) = ψ 0 + L f g T ( T - T ) T < T ψ 0 T T ,

where ψ0 is the matric potential due to saturation, Lf= 3.4 × 105J kg−1 is the specific latent heat of fusion of water, and g= 9.80665 m s−2 is downward acceleration due to gravity. The critical temperature T is expressed as a function of the melting temperature Tm and the matric potential:

(B2) T = T m + g T m L f ψ 0 .

The volumetric liquid water content then follows from the well-known van Genuchten equation (van Genuchten1980):

(B3) θ w = θ r + ( θ s - θ r ) ( 1 + ( - α ψ ( T ) ) n ) - m ,

where α(0,) [m−1], n[1,), and m=1-1n are empirical fitting parameters of the van Genuchten curve which depend on soil type, and θr and θs are the residual and saturated water contents, respectively. We use the log-transform of the strictly positive van Genuchten parameters in the inversion procedure to (i) allow for the use of Gaussian priors in unconstrained space and (ii) to improve the numerical properties of the likelihood which, as should be apparent from Eq. (B3), exhibit strong nonlinearities in the derivatives with respect to α and n.

B2.2 Numerical scheme

The spatial gradients in Eq. (3) are approximated using the method of lines:

(B4) t H i ( z , t ) = 1 Δ i [ k i ( T i + 1 - T i ) δ i - k i - 1 ( T i - T i - 1 ) δ i - 1 ] ,

where Δi is the length of the ith inner grid cell, and δi is the distance between the midpoints of cells i and i+1. Since soil properties and saturation conditions are held constant in our simulations, we reduce the computational requirements of the model by pre-computing the volumetric unfrozen water content as a function of enthalpy, θw(H), for subzero temperatures in the range 50 to 0 C with interpolation knots chosen such that the maximum permissible error of the approximation is controlled to be less than 0.001 % or 0.1 % of the volumetric fraction.

Table B5Constituent thermal property values used for all stratigraphy layers and all sites.

Download Print Version | Download XLSX

The bulk, temperature-dependent thermal conductivity k(T) is parameterized following (Cosenza et al.2003) as the inverse quadratic mean of the constituents:

(B5) k ( T ) = ( k w θ w ( T ) + k i θ i ( T ) + k a θ a + k m θ m + k o θ o ) 2 ,

where k [Wm-1K-1] refers to the thermal conductivity of the constituent material. The bulk heat capacity is similarly computed as a simple weighted average:

(B6) C ( T ) = 1 5 ( c w θ w ( T ) + c i θ i ( T ) + c a θ a + c m θ m + c o θ o ) ,

where c values are the constituent heat capacities. Since we already vary the constituent fractions θ in our inversion procedure, we hold these values for the constituent thermal properties constant in all of our simulations. The values are shown in Table B5.

For time integration, we use an implementation of the second-order, explicit, strong-stability-preserving time stepping scheme of Shu and Osher (1988) provided by the Julia package, OrdinaryDiffEq.jl (Rackauckas and Nie2017). The time step size is set according to a suitable Courant–Friedrichs–Lewy (CFL) condition (Courant et al.1967) for the nonlinear heat equation:

(B7) d t max = C min i H T i k i - 1 Δ i 2 ,

where HTi [JK-1m-3] is the apparent heat capacity at grid cell i, ki [Wm-1K-1] is the thermal conductivity, Δi [m] is the vertical length of the grid cell, and C=12 is the dimensionless Courant number which scales the step size to maintain stability.

B3 Inversion method

B3.1 Posterior predictions

From Eq. (7), the predictive distribution over the model states can be formulated as

(B8) p ( x 1 : τ | y 1 : τ , s 1 : τ ) = ϕ Φ p ( x 1 : τ | ϕ , s 1 : τ ) × p ( ϕ | y 1 : τ , s 1 : τ ) d ϕ ,

with the process model p(x1:τ|ϕ,s1:τ)=δ(x1:τ=f(ϕ,s1:τ)), where f is the transient heat conduction model (described in Sect. 3.3) parameterized by ϕ, and δ is the Dirac delta distribution. Note that this formulation of the problem departs somewhat from typical data assimilation procedures which tend to focus more on probabilistic modeling of the state x, with parametric uncertainty in ϕ often neglected (Geer2021; Caers2018). In contrast, our approach amounts to a kind of “perfect model” assumption where the dynamics are assumed to be error-free. This allows us to avoid the need to develop a stochastic description of p(x1:τ|ϕ,s1:τ), while still being able to produce a distribution over plausible unobserved states via the parametric uncertainty induced by the prior. In principle, model error can be implicitly accounted for in the likelihood p(y1:τ|x1:τ,ϕ,s1:τ). However, as discussed in Sect. 5.5, EKS assumes the noise scale of the likelihood to be fixed and known a priori, which is generally insufficient for nontrivial inverse problems. As a result, the predictive distribution described by Eq. (B8) cannot be guaranteed to be well calibrated.

Estimating the predictive distribution in Eq. (B8) requires integration over the parameter space Φ weighted by the posterior distribution p(ϕ|y1:τ,s1:τ) which is generally intractable. A Monte Carlo approximation to the integral can, however, be obtained by using sampling algorithms such as EKS to obtain samples from the posterior distribution.

B3.2 EKS problem setup

Since we aim to constrain the parameters ϕ of the heat conduction model f based on observed borehole temperatures, we choose gT in Eq. (8) to be a simple aggregation function which computes mean annual ground temperatures (MAGTs) from the daily temperature profiles produced by f. For each site, we compute the modeled MAGT at each depth where a borehole sensor is present, linearly interpolating between grid cells where necessary. The resulting predicted temperatures form a matrix Tpred=gTf(ϕ)RN×d, where N is the number of years and d the number of borehole sensor depths at which data are available. Model outputs at all points in time and space where no data are present are thus not explicitly included in the likelihood calculation, though it is important to note that they are still included implicitly through conditional dependence of the evaluation points on the time dynamics of Eq. 3. We use the observed mean annual ground temperature from each borehole site as the observations, i.e., Tobs in Eq. (8), in the EKS algorithm. We specify the observation covariance ΣT as a diagonal matrix with the posterior predictive variances obtained from our trend analysis (see Sects. 4.1 and B1) along the diagonal and a minimum variance of 0.01 K2 (based on the temperature resolution of the thermistor chain, 0.1 K, according to manufacturer specifications). The diagonal structure of ΣT assumes the observation noise to be both independent and uncorrelated across time and space.

The prior distribution over model parameters p(ϕ) is of crucial importance to any Bayesian inference procedure. Some parameters in the heat transfer model, such as soil composition, will cause the resulting optimization problem on ϕ to be under-constrained, since there may be more than one possible combination of soil components which have similar thermal properties. Additionally, incorporating prior knowledge about plausible parameter values allows us to reduce the amount of computational effort wasted on physically implausible or incoherent model configurations that may arise from random sampling. We define the priors for the heat transfer model parameters ϕ in their native, constrained parameter space Φ in order to more easily incorporate physically meaningful information about each site. More details on choices of priors for each site are in Appendix B3.3. We construct a diagonal, multivariate Gaussian approximation of our priors in the unconstrained parameter space via Laplace approximation of the density induced by γ(ψ). Note that this introduces some error into the EKS representation of the posterior, since the resulting Gaussian distribution is only an approximation of the transformed density.

B3.3 Selection of prior distributions

There are broadly three categories of parameters which are considered in our analysis: soil parameters (see Appendix B2.1), n factor parameters, and initial condition parameters. Separate soil parameters are defined for each stratigraphy layer. Soil freeze curve parameters (melting temperature Tm and van Genuchten parameters log α, log n) as well as porosity are varied in all layers, while organic content is only varied in the uppermost soil layer and is generally assumed to be zero or negligible in all other layers. Ground ice content is also varied only in layers where it is suspected to be present. Saturation of pore-space water and ice is varied as a single, global parameter used by each layer, excluding the bottom-most layer, which is held constant; since saturation is temporally constant throughout the simulation, this can be seen as effectively a bulk or average saturation level for the soil column.

We based the stratigraphy priors for Samoylov, Bayelva, and Barrow on previous modeling studies, Nitzbon et al. (2021), Zweigel et al. (2021), and Romanovsky and Osterkamp (2000), respectively. For Parson's Lake, we designed the stratigraphy based on the soil core analysis for the borehole site presented by Wolfe et al. (2010). We chose the priors for the van Genuchten parameters by fitting empirical freeze curves to temperature and volumetric water content measurements available for both the Samoylov Island and Bayelva sites (Boike et al.2018a, 2019). The mean and standard deviation of the log-transformed parameters are set according to the fitted posterior statistics, with the variance increased by a factor of 100 to roughly account for uncertainty in extrapolation. We omit the details of this fitting procedure here for brevity, but the relevant code and results are included in the supplementary material. For Barrow, we used the same van Genuchten parameter priors derived from the 38 cm water content sensor at the Bayelva, since the freeze curve visually agreed with the one presented in Romanovsky and Osterkamp (2000), and both sites are known to have silty soil in the upper half meter (Boike et al.2018a). We set these priors for Parson's Lake based on the reported soil types in each layer (Wolfe et al.2010) using the Samoylov priors in the upper sandy layers and the Bayelva silty clay priors for the deeper layers. We stress that, in all cases, there is currently no systematic procedure used to derive the prior distributions; the process is largely based on subjective evaluation and iterative validation. Furthermore, the lack of stratigraphy information for the deeper layers is a substantial source of uncertainty in our analysis. While the variability induced by sampling from the prior distributions accounts for some of this uncertainty, there is a lot of room for improvement in the development of more rigorous and principled methods to derive prior distributions for soil stratigraphies and their corresponding parameterizations.

Code and data availability

The preprocessed/collated data files (excluding those from Parson's Lake), in addition to the code for this study, are accessible online at (Groenke et al.2023). The code for the transient heat conduction model is also available online at (Groenke et al.2022). The borehole data from the Parson's Lake borehole can be requested directly from the Geological Survey of Canada, Natural Resources Canada (Smith et al.2018; Wolfe et al.2010). The original borehole data from the Samoylov Island and Bayelva sites can be retrieved from the open-access Pangaea data store (Boike et al.2018b, 2017). The original borehole data from the Barrow site can be retrieved from the University of Alaska Fairbanks Permafrost Laboratory website (Bar2021).


The supplement related to this article is available online at:

Author contributions

BG, ML, and SW designed the study and modeling scheme; BG designed and authored the model and analysis code, carried out all numerical experiments, and wrote the paper with support from all co-authors. JB and GG supervised the project. All authors discussed the results and contributed to the final version of the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We are grateful to Caroline Duchense and Sharon Smith (Natural Resources Canada) for providing us with the borehole temperature record from Parson's Lake and to Adam Theisen (Argonne National Laboratory) for providing climatology data from the Northern Slope of Alaska. We are additionally grateful to both Thomas Hermans and the anonymous reviewer of this work for their constructive criticism and feedback. We would also like to thank Inge Grünberg (Alfred Wegener Institute) for her valuable consultation and feedback.

Financial support

This research has been supported by the Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, and the Helmholtz Einstein International Berlin Research School in Data Science (grant no. HIDSS-0001).

The article processing charges for this open-access publication were covered by the Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research (AWI).

Review statement

This paper was edited by Adrian Flores Orozco and reviewed by Thomas Hermans and one anonymous referee.


Allen, D. M., Michel, F. A., and Judge, A. S.: The Permafrost Regime in the Mackenzie Delta, Beaufort Sea Region, N. W. T. and Its Significance to the Reconstruction of the Palaeoclimatic History, J. Quaternary Sci., 3, 3–13,, 1988. a

Barrow 2: N. Meadow Lake No.2/NML-2, (last access: 6 July 2022), 2021. a, b, c

Berliner, L. M.: Physical-Statistical Modeling in Geophysics, J. Geophys. Res.-Atmos., 108, 8776,, 2003. a, b

Bishop, C.: Pattern Recognition and Machine Learning, 1st Edn., vol. 4 of Information Science and Statistics, Springer-Verlag, New York, ISBN: 978-0-387-31073-2, 2006. a

Biskaborn, B. K., Lanckman, J.-P., Lantuit, H., Elger, K., Streletskiy, D. A., Cable, W. L., and Romanovsky, V. E.: The new database of the Global Terrestrial Network for Permafrost (GTN-P), Earth Syst. Sci. Data, 7, 245–259,, 2015. a, b

Biskaborn, B. K., Smith, S. L., Noetzli, J., Matthes, H., Vieira, G., Streletskiy, D. A., Schoeneich, P., Romanovsky, V. E., Lewkowicz, A. G., Abramov, A., Allard, M., Boike, J., Cable, W. L., Christiansen, H. H., Delaloye, R., Diekmann, B., Drozdov, D., Etzelmüller, B., Grosse, G., Guglielmin, M., Ingeman-Nielsen, T., Isaksen, K., Ishikawa, M., Johansson, M., Johannsson, H., Joo, A., Kaverin, D., Kholodov, A., Konstantinov, P., Kröger, T., Lambiel, C., Lanckman, J.-P., Luo, D., Malkova, G., Meiklejohn, I., Moskalenko, N., Oliva, M., Phillips, M., Ramos, M., Sannel, A. B. K., Sergeev, D., Seybold, C., Skryabin, P., Vasiliev, A., Wu, Q., Yoshikawa, K., Zheleznyak, M., and Lantuit, H.: Permafrost Is Warming at a Global Scale, Nat. Commun., 10, 264,, 2019. a, b, c, d, e, f

Boike, J., Juszak, I., Lange, S., Chadburn, S., Burke, E. J., Overduin, P. P., Roth, K., Ippisch, O., Bornemann, N., Stern, L., Gouttevin, I., Hauber, E., and Westermann, S.: Measurements in Soil and Air at Bayelva Station, PANGAEA [data set],, 2017. a

Boike, J., Juszak, I., Lange, S., Chadburn, S., Burke, E., Overduin, P. P., Roth, K., Ippisch, O., Bornemann, N., Stern, L., Gouttevin, I., Hauber, E., and Westermann, S.: A 20-year record (1998–2017) of permafrost, active layer and meteorological conditions at a high Arctic permafrost research site (Bayelva, Spitsbergen), Earth Syst. Sci. Data, 10, 355–390,, 2018a. a, b, c, d, e, f, g

Boike, J., Nitzbon, J., Anders, K., Grigoriev, M. N., Bolshiyanov, D. Y., Langer, M., Lange, S., Bornemann, N., Morgenstern, A., Schreiber, P., Wille, C., Chadburn, S., Gouttevin, I., and Kutzbach, L.: Measurements in Soil and Air at Samoylov Station (2002–2018), PANGAEA [data set],, 2018b. a

Boike, J., Nitzbon, J., Anders, K., Grigoriev, M., Bolshiyanov, D., Langer, M., Lange, S., Bornemann, N., Morgenstern, A., Schreiber, P., Wille, C., Chadburn, S., Gouttevin, I., Burke, E., and Kutzbach, L.: A 16-year record (2002–2017) of permafrost, active-layer, and meteorological conditions at the Samoylov Island Arctic permafrost research site, Lena River delta, northern Siberia: an opportunity to validate remote-sensing data and land surface, snow, and permafrost models, Earth Syst. Sci. Data, 11, 261–299,, 2019. a, b, c, d, e, f, g, h

Brewer, M. C.: Some Results of Geothermal Investigations of Permafrost in Northern Alaska, Eos T. Am. Geophys. Un., 39, 19–26,, 1958. a

Brown, J., Ferrians Jr, O. J., Heginbottom, J. A., and Melnikov, E. S.: Circum-Arctic Map of Permafrost and Ground-Ice Conditions, US Geological Survey Reston, 1997. a

Burgess, M. M., Smith, S. L., Brown, J., Romanovsky, V., and Hinkel, K.: Global Terrestrial Network for Permafrost (GTNet-P): Permafrost Monitoring Contributing to Global Climate Observations, Tech. Rep. 2000-E14,, 2000. a

Burn, C. R. and Kokelj, S. V.: The Environment and Permafrost of the Mackenzie Delta Area, Permafrost Periglac., 20, 83–105,, 2009. a, b

Caers, J.: Bayesianism in the Geosciences, in: Handbook of Mathematical Geosciences: Fifty Years of IAMG, edited by: Daya Sagar, B., Cheng, Q., and Agterberg, F., Springer International Publishing, Cham, 527–566,, 2018. a, b, c

Chen, L., Aalto, J., and Luoto, M.: Significant Shallow–Depth Soil Warming over Russia during the Past 40 Years, Global Planet. Change, 197, 103394,, 2021. a, b, c

Cleary, E., Garbuno-Inigo, A., Lan, S., Schneider, T., and Stuart, A. M.: Calibrate, Emulate, Sample, J. Comput. Physics, 424, 109716,, 2021. a, b

Constantinou, N. C., Lopez-Gomez, I., Iñigo, A. G., Bolewski, J., odunbar, Howland, M., Wu, J., Jackson, T., (Zhengyu)Huang, D., Jaruga, A., Ervik, H. L. L., Monticone, P., Chor, T., and bielim: CliMA/EnsembleKalmanProcesses.Jl: V0.8.0, Zenodo [code],, 2022. a

Coon, E., Svyatsky, D., Jan, A., Kikinzon, E., Berndt, M., Atchley, A., Harp, D., Manzini, G., Shelef, E., Lipnikov, K., Garimella, R., Xu, C., Moulton, D., Karra, S., and Painter, S.: Advanced Terrestrial Simulator, USDOE Office of Science (SC), Biological and Environmental Research (BER) (SC-23),, 2019. a

Cosenza, P., Guérin, R., and Tabbagh, A.: Relationship between Thermal Conductivity and Water Content of Soils Using Numerical Modelling, Eur. J. Soil Sci., 54, 581–588,, 2003. a

Courant, R., Friedrichs, K., and Lewy, H.: On the Partial Difference Equations of Mathematical Physics, IBM J. Res. Dev., 11, 215–234, 1967. a

Cranmer, K., Brehmer, J., and Louppe, G.: The Frontier of Simulation-Based Inference, P. Natl. Acad. Sci. USA, 117, 30055–30062,, 2020. a, b

Cuesta-Valero, F. J., Beltrami, H., García-García, A., Krinner, G., Langer, M., MacDougall, A. H., Nitzbon, J., Peng, J., von Schuckmann, K., Seneviratne, S. I., Thiery, W., Vanderkelen, I., and Wu, T.: Continental heat storage: contributions from the ground, inland waters, and permafrost thawing, Earth Syst. Dynam., 14, 609–627,, 2023. a

Dall'Amico, M., Endrizzi, S., Gruber, S., and Rigon, R.: A robust and energy-conserving model of freezing variably-saturated soil, The Cryosphere, 5, 469–484,, 2011. a, b

Doucet, A., de Freitas, N., and Gordon, N.: An Introduction to Sequential Monte Carlo Methods, in: Sequential Monte Carlo Methods in Practice, edited by Doucet, A., de Freitas, N., and Gordon, N., Statistics for Engineering and Information Science, Springer, New York, NY,, 3–14, 2001. a

Duane, S. and Kogut, J. B.: The Theory of Hybrid Stochastic Algorithms, Nucl. Phys. B, 275, 398–420,, 1986. a

ECCC: Historical Data – Climate – Environment and Climate Change Canada, (last access: 8 July 2022), 2011. a

Evensen, G., Vossepoel, F. C., and van Leeuwen, P. J.: Data Assimilation Fundamentals: A Unified Formulation of the State and Parameter Estimation Problem, Springer Nature,, 2022. a

Franks, S. W. and Beven, K. J.: Bayesian Estimation of Uncertainty in Land Surface-Atmosphere Flux Predictions, J. Geophys. Res.-Atmos., 102, 23991–23999,, 1997. a

Garbuno-Inigo, A., Hoffmann, F., Li, W., and Stuart, A. M.: Interacting Langevin Diffusions: Gradient Structure and Ensemble Kalman Sampler, SIAM J. Appl. Dyn. Syst., 19, 412–441,, 2020. a, b, c, d, e

Garnello, A., Marchenko, S., Nicolsky, D., Romanovsky, V., Ledman, J., Celis, G., Schädel, C., Luo, Y., and Schuur, E. a. G.: Projecting Permafrost Thaw of Sub-Arctic Tundra With a Thermodynamic Model Calibrated to Site Measurements, J. Geophys. Res.-Bioge., 126, e2020JG006218,, 2021. a, b

Geer, A. J.: Learning Earth System Models from Observations: Machine Learning or Data Assimilation?, Philos. T. R. Soc. A, 379, 20200089,, 2021. a

Gopalan, G., Hrafnkelsson, B., Aðalgeirsdóttir, G., Jarosch, A. H., and Pálsson, F.: A Bayesian hierarchical model for glacial dynamics based on the shallow ice approximation and its evaluation using analytical solutions, The Cryosphere, 12, 2229–2248,, 2018. a

Gregory, W., Lawrence, I. R., and Tsamados, M.: A Bayesian approach towards daily pan-Arctic sea ice freeboard estimates from combined CryoSat-2 and Sentinel-3 satellite observations, The Cryosphere, 15, 2857–2871,, 2021. a

Grigoriev, N. F.: The Temperature of Permafrost in the Lena Delta Basin–Deposit Conditions and Properties of the Permafrost in Yakutia, Yakutsk, 2, 97–101, 1960. a

Groenke, B., Nitzbon, J., and Langer, M.: CryoGrid.Jl: V0.10.3, Zenodo [code],, 2022. a

Groenke, B., Langer, M., Gallego, G., and Boike, J.: Investigating the Thermal State of Permafrost with Bayesian Inverse Modeling of Heat Transfer – Data and Code Repository, Zenodo [code and data set],, 2023. a

Harp, D. R., Atchley, A. L., Painter, S. L., Coon, E. T., Wilson, C. J., Romanovsky, V. E., and Rowland, J. C.: Effect of soil property uncertainties on permafrost thaw projections: a calibration-constrained analysis, The Cryosphere, 10, 341–358,, 2016. a

Hermans, T., Nguyen, F., Klepikova, M., Dassargues, A., and Caers, J.: Uncertainty Quantification of Medium-Term Heat Storage From Short-Term Geophysical Experiments Using Bayesian Evidential Learning, Water Resour. Res., 54, 2931–2948,, 2018. a, b, c

Hermoso de Mendoza, I., Beltrami, H., MacDougall, A. H., and Mareschal, J.-C.: Lower boundary conditions in land surface models – effects on the permafrost and the carbon pools: a case study with CLM4.5, Geosci. Model Dev., 13, 1663–1683,, 2020. a

Hinzman, L. D., Bettez, N. D., Bolton, W. R., Chapin, F. S., Dyurgerov, M. B., Fastie, C. L., Griffith, B., Hollister, R. D., Hope, A., Huntington, H. P., Jensen, A. M., Jia, G. J., Jorgenson, T., Kane, D. L., Klein, D. R., Kofinas, G., Lynch, A. H., Lloyd, A. H., McGuire, A. D., Nelson, F. E., Oechel, W. C., Osterkamp, T. E., Racine, C. H., Romanovsky, V. E., Stone, R. S., Stow, D. A., Sturm, M., Tweedie, C. E., Vourlitis, G. L., Walker, M. D., Walker, D. A., Webber, P. J., Welker, J. M., Winker, K. S., and Yoshikawa, K.: Evidence and Implications of Recent Climate Change in Northern Alaska and Other Arctic Regions, Climatic Change, 72, 251–298,, 2005. a

Hoffman, M. D. and Gelman, A.: The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo, J. Mach Learn. Res., 15, 31,, 2014. a, b

Humlum, O.: Holocene Permafrost Aggradation in Svalbard, Geological Society, London, Special Publications, 242, 119–129,, 2005. a

Iglesias, M. A., Law, K. J. H., and Stuart, A. M.: Ensemble Kalman Methods for Inverse Problems, Inverse Probl., 29, 045001,, 2013. a

Isaksen, K., Nordli, Ø., Ivanov, B., Køltzow, M. A. Ø., Aaboe, S., Gjelten, H. M., Mezghani, A., Eastwood, S., Førland, E., Benestad, R. E., Hanssen-Bauer, I., Brækkan, R., Sviashchennikov, P., Demin, V., Revina, A., and Karandasheva, T.: Exceptional Warming over the Barents Area, Sci. Rep.-UK, 12, 9371,, 2022. a

Jafarov, E. E., Marchenko, S. S., and Romanovsky, V. E.: Numerical modeling of permafrost dynamics in Alaska using a high spatial resolution dataset, The Cryosphere, 6, 613–624,, 2012. a

Jorgenson, M. T., Shur, Y. L., and Pullman, E. R.: Abrupt Increase in Permafrost Degradation in Arctic Alaska, Geophys. Res. Lett., 33, L02503,, 2006. a

Kantas, N., Beskos, A., and Jasra, A.: Sequential Monte Carlo Methods for High-Dimensional Inverse Problems: A Case Study for the Navier–Stokes Equations, SIAM/ASA Journal on Uncertainty Quantification, 2, 464–489,, 2014. a

Kneier, F., Overduin, P. P., Langer, M., Boike, J., and Grigoriev, M. N.: Borehole Temperature Reconstructions Reveal Differences in Past Surface Temperature Trends for the Permafrost in the Laptev Sea Region, Russian Arctic, arktos, 4, 1–17,, 2018. a

Koopmans, R. W. R. and Miller, R. D.: Soil Freezing and Soil Water Characteristic Curves, Soil Sci. Soc. Am. J., 30, 680–685, 1966. a

Lange, K. L., Little, R. J., and Taylor, J. M.: Robust Statistical Modeling Using the t Distribution, J. Am. Stat. Assoc., 84, 881–896, 1989. a

Langer, M., Westermann, S., Muster, S., Piel, K., and Boike, J.: The surface energy balance of a polygonal tundra site in northern Siberia – Part 1: Spring to fall, The Cryosphere, 5, 151–171,, 2011a. a

Langer, M., Westermann, S., Muster, S., Piel, K., and Boike, J.: The surface energy balance of a polygonal tundra site in northern Siberia – Part 2: Winter, The Cryosphere, 5, 509–524,, 2011b. a

Langer, M., Nitzbon, J., Groenke, B., Assmann, L.-M., Schneider von Deimling, T., Stuenzi, S. M., and Westermann, S.: The evolution of Arctic permafrost over the last three centuries, EGUsphere [preprint],, 2022. a, b

Lunardini, V. J.: A Correlation of N-Factors, in: Applied Techniques for Cold Environments, ASCE, 233–244, 1978. a

Moradkhani, H., Hsu, K.-L., Gupta, H., and Sorooshian, S.: Uncertainty Assessment of Hydrologic Model States and Parameters: Sequential Data Assimilation Using the Particle Filter, Water Resour. Res., 41, 5,, 2005. a

Muñoz-Sabater, J., Dutra, E., Agustí-Panareda, A., Albergel, C., Arduini, G., Balsamo, G., Boussetta, S., Choulga, M., Harrigan, S., Hersbach, H., Martens, B., Miralles, D. G., Piles, M., Rodríguez-Fernández, N. J., Zsoter, E., Buontempo, C., and Thépaut, J.-N.: ERA5-Land: a state-of-the-art global reanalysis dataset for land applications, Earth Syst. Sci. Data, 13, 4349–4383,, 2021. a

Nelson, F. E., Anisimov, O. A., and Shiklomanov, N. I.: Climate Change and Hazard Zonation in the Circum-Arctic Permafrost Regions, Nat. Hazards, 26, 203–225,, 2002. a

Nelson, F. E., Shiklomanov, N. I., Streletskiy, D. A., Romanovsky, V. E., Yoshikawa, K., Hinkel, K. M., and Brown, J.: A Permafrost Observatory at Barrow, Alaska: Longterm Observations of Active-Layer Thickness and Permafrost Temperature, in: Ninth International Conference on Permafrost, 2, 1267–1280, 2008. a

Nicolsky, D. J. and Romanovsky, V. E.: Modeling Long-Term Permafrost Degradation, J. Geophys. Res.-Earth, 123, 1756–1771,, 2018. a, b, c, d, e

Nicolsky, D. J., Romanovsky, V. E., Alexeev, V. A., and Lawrence, D. M.: Improved Modeling of Permafrost Dynamics in a GCM Land-Surface Scheme: Modeling Permafrost Dynamics in Alaska by CLM3, Geophys. Res. Lett., 34, L08501,, 2007. a

Nicolsky, D. J., Romanovsky, V. E., and Panteleev, G. G.: Estimation of Soil Thermal Properties Using In-Situ Temperature Measurements in the Active Layer and Permafrost, Cold Reg. Sci. Technol., 55, 120–129,, 2009. a

Nitzbon, J., Langer, M., Martin, L. C. P., Westermann, S., Schneider von Deimling, T., and Boike, J.: Effects of multi-scale heterogeneity on the simulated evolution of ice-rich permafrost lowlands under a warming climate, The Cryosphere, 15, 1399–1422,, 2021. a

Nitzbon, J., Krinner, G., Schneider Von Deimling, T., Werner, M., and Langer, M.: First Quantification of the Permafrost Heat Sink in the Earth's Climate System, Geophys. Res. Lett., 50, e2022GL102053,, 2023. a, b

Noh, S. J., Tachikawa, Y., Shiiba, M., and Kim, S.: Applying sequential Monte Carlo methods into a distributed hydrologic model: lagged particle filtering approach with regularization, Hydrol. Earth Syst. Sci., 15, 3237–3251,, 2011. a

Obu, J., Westermann, S., Bartsch, A., Berdnikov, N., Christiansen, H. H., Dashtseren, A., Delaloye, R., Elberling, B., Etzelmüller, B., Kholodov, A., Khomutov, A., Kääb, A., Leibman, M. O., Lewkowicz, A. G., Panda, S. K., Romanovsky, V., Way, R. G., Westergaard-Nielsen, A., Wu, T., Yamkhin, J., and Zou, D.: Northern Hemisphere Permafrost Map Based on TTOP Modelling for 2000–2016 at 1 km2 Scale, Earth-Sci. Rev., 193, 299–316,, 2019. a

Park, H., Fedorov, A. N., Zheleznyak, M. N., Konstantinov, P. Y., and Walsh, J. E.: Effect of Snow Cover on Pan-Arctic Permafrost Thermal Regimes, Clim. Dynam., 44, 2873–2895,, 2015. a

Piani, C., Weedon, G. P., Best, M., Gomes, S. M., Viterbo, P., Hagemann, S., and Haerter, J. O.: Statistical Bias Correction of Global Simulated Daily Precipitation and Temperature for the Application of Hydrological Models, J. Hydrol., 395, 199–215,, 2010. a

Qu, Y., Wang, J., Wan, H., Li, X., and Zhou, G.: A Bayesian Network Algorithm for Retrieving the Characterization of Land Surface Vegetation, Remote Sens. Environ., 112, 613–622,, 2008. a

Rackauckas, C. and Nie, Q.: DifferentialEquations.Jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia, Journal of Open Research Software, 5, 15,, 2017. a

Riseborough, D., Shiklomanov, N., Etzelmüller, B., Gruber, S., and Marchenko, S.: Recent advances in permafrost modelling, Permafrost Periglac., 19, 137–156,, 2008. a, b

Riseborough, D. W.: Soil Latent Heat as a Filter of the Climate Signal in Permafrost, in: Proceedings of the Fifth Canadian Permafrost Conference, Collection Nordicana, vol. 54, pp. 199–205, Citeseer, 1990. a, b, c, d

Romanovsky, V., Burgess, M., Smith, S., Yoshikawa, K., and Brown, J.: Permafrost Temperature Records: Indicators of Climate Change, Eos T. Am. Geophys. Un., 83, 589–594,, 2002. a

Romanovsky, V. E. and Osterkamp, T. E.: Thawing of the Active Layer on the Coastal Plain of the Alaskan Arctic, Permafrost Periglac., 8, 1–22,<1::AID-PPP243>3.0.CO;2-U, 1997. a

Romanovsky, V. E. and Osterkamp, T. E.: Effects of Unfrozen Water on Heat and Mass Transport Processes in the Active Layer and Permafrost, Permafrost Periglac., 11, 219–239,<219::AID-PPP352>3.0.CO;2-7, 2000. a, b, c, d, e

Romanovsky, V. E., Sazonova, T. S., Balobaev, V. T., Shender, N. I., and Sergueev, D. O.: Past and Recent Changes in Air and Permafrost Temperatures in Eastern Siberia, Global Planet. Change, 56, 399–413,, 2007. a, b

Romanovsky, V. E., Smith, S. L., and Christiansen, H. H.: Permafrost Thermal State in the Polar Northern Hemisphere during the International Polar Year 2007–2009: A Synthesis, Permafrost Periglac., 21, 106–116,, 2010. a

Roth, K. and Boike, J.: Quantifying the Thermal Dynamics of a Permafrost Site near Ny-Ålesund, Svalbard, Water Resour. Res., 37, 2901–2914,, 2001. a

Rubin, D. B.: Bayesianly Justifiable and Relevant Frequency Calculations for the Applied Statistician, Ann. Stat., 12, 1151–1172, 1984. a

Schneider von Deimling, T., Lee, H., Ingeman-Nielsen, T., Westermann, S., Romanovsky, V., Lamoureux, S., Walker, D. A., Chadburn, S., Trochim, E., Cai, L., Nitzbon, J., Jacobi, S., and Langer, M.: Consequences of permafrost degradation for Arctic infrastructure – bridging the model gap between regional and engineering scales, The Cryosphere, 15, 2451–2471,, 2021. a

Schuur, E. A. and Mack, M. C.: Ecological Response to Permafrost Thaw and Consequences for Local and Global Ecosystem Services, Annu. Rev. Ecol. Evol. S., 49, 279–301,, 2018. a

Serreze, M. C. and Francis, J. A.: The Arctic Amplification Debate, Climatic Change, 76, 241–264,, 2006. a

Shiklomanov, N. I., Streletskiy, D. A., Nelson, F. E., Hollister, R. D., Romanovsky, V. E., Tweedie, C. E., Bockheim, J. G., and Brown, J.: Decadal Variations of Active-Layer Thickness in Moisture-Controlled Landscapes, Barrow, Alaska, J. Geophys. Res.-Biogeo., 115, G00I04,, 2010. a, b

Shu, C.-W. and Osher, S.: Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing Schemes, J. Comput. Phys., 77, 439–471, 1988. a

Sisson, S. A., Fan, Y., and Tanaka, M. M.: Sequential Monte Carlo without Likelihoods, P. Natl. Acad. Sci. USA, 104, 1760–1765,, 2007. a

Sisson, S. A., Fan, Y., and Beaumont, M.: Handbook of Approximate Bayesian Computation, CRC Press,, 2018. a

Smith, S., Romanovsky, V., Lewkowicz, A., Burn, C., Allard, M., Clow, G., Yoshikawa, K., and Throop, J.: Thermal State of Permafrost in North America: A Contribution to the International Polar Year, Permafrost Periglac., 21, 117–135,, 2010. a

Smith, S. L., O'Neill, H. B., Isaksen, K., Noetzli, J., and Romanovsky, V. E.: The Changing Thermal State of Permafrost, Nature Reviews Earth & Environment, 3, 10–23,, 2022. a, b, c, d

Smith, S. L., Chartrand, J., and Duchesne, C.: Report on 2017 Field Activities and Collection of Ground-Thermal and Active-Layer Data in the Mackenzie Corridor, Northwest Territories, Geological Survey of Canada Open File 8492, Geological Survey of Canada, Ottawa, ON, 109 pp, 2018. a, b, c

Theisen, A.: ARM-Climatologies North Slope Alaska, Zenodo,, 2022. a

Thibaut, R., Compaire, N., Lesparre, N., Ramgraber, M., Laloy, E., and Hermans, T.: Comparing Well and Geophysical Data for Temperature Monitoring within a Bayesian Experimental Design Framework, Water Resour. Res., 58, e2022WR033045,, 2022.  a, b

Tubini, N., Gruber, S., and Rigon, R.: A method for solving heat transfer with phase change in ice or soil that allows for large time steps while guaranteeing energy conservation, The Cryosphere, 15, 2541–2568,, 2021. a, b

van Genuchten, M. T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., 44, 892–898,, 1980. a, b, c

Verjans, V., Leeson, A. A., Nemeth, C., Stevens, C. M., Kuipers Munneke, P., Noël, B., and van Wessem, J. M.: Bayesian calibration of firn densification models, The Cryosphere, 14, 3017–3032,, 2020. a

Wainwright, H. M., Liljedahl, A. K., Dafflon, B., Ulrich, C., Peterson, J. E., Gusmeroli, A., and Hubbard, S. S.: Mapping snow depth within a tundra ecosystem using multiscale observations and Bayesian methods, The Cryosphere, 11, 857–875,, 2017. a

Westermann, S., Boike, J., Langer, M., Schuler, T. V., and Etzelmüller, B.: Modeling the impact of wintertime rain events on the thermal regime of permafrost, The Cryosphere, 5, 945–959,, 2011. a

Westermann, S., Langer, M., Boike, J., Heikenfeld, M., Peter, M., Etzelmüller, B., and Krinner, G.: Simulating the thermal regime and thaw processes of ice-rich permafrost ground with the land-surface model CryoGrid 3, Geosci. Model Dev., 9, 523–546,, 2016. a

Westermann, S., Ingeman-Nielsen, T., Scheer, J., Aalstad, K., Aga, J., Chaudhary, N., Etzelmüller, B., Filhol, S., Kääb, A., Renette, C., Schmidt, L. S., Schuler, T. V., Zweigel, R. B., Martin, L., Morard, S., Ben-Asher, M., Angelopoulos, M., Boike, J., Groenke, B., Miesner, F., Nitzbon, J., Overduin, P., Stuenzi, S. M., and Langer, M.: The CryoGrid community model (version 1.0) – a multi-physics toolbox for climate-driven simulations in the terrestrial cryosphere, Geosci. Model Dev., 16, 2607–2647,, 2023. a, b, c

Wolfe, S. A., Smith, S. L., Chartrand, J., Kokelj, S., Palmer, M., and Stevens, C. W.: Geotechnical Database and Descriptions of Permafrost Monitoring Sites Established 2006-2010 in the Northern Mackenzie Corridor, Northwest Territories, Tech. Rep. 6677,, 2010. a, b, c, d

Zweigel, R. B., Westermann, S., Nitzbon, J., Langer, M., Boike, J., Etzelmüller, B., and Vikhamar Schuler, T.: Simulating Snow Redistribution and Its Effect on Ground Surface Temperature at a High-Arctic Site on Svalbard, J. Geophys, Res.-Earth, 126, e2020JF005673,, 2021. a

Short summary
It is now well known from long-term temperature measurements that Arctic permafrost, i.e., ground that remains continuously frozen for at least 2 years, is warming in response to climate change. Temperature, however, only tells half of the story. In this study, we use computer modeling to better understand how the thawing and freezing of water in the ground affects the way permafrost responds to climate change and what temperature trends can and cannot tell us about how permafrost is changing.