PERICLIMv1.0: a model deriving palaeo-air temperatures from thaw depth in past permafrost regions

Periglacial features, such as various kinds of patterned ground, cryoturbations, frost wedges, solifluction structures, and blockfields, are among the most common relics of cold climate periods, which repetitively occurred throughout the Quaternary. As such, they are widespread archives of past environmental conditions. Climate controls on the development of most periglacial features, however, remain poorly known, and thus empirical palaeo-climate reconstructions based on them have limited validity. This study presents and evaluates a simple new inverse modelling scheme called PERICLIMv1.0 (PERIglacial CLIMate) that derives palaeo-air temperature characteristics related to the palaeo-active-layer thickness, which can be recognized using many relict periglacial features found in past permafrost regions. The evaluation against modern temperature records showed that the model reproduces air temperature characteristics with average errors ≤1.3 ^∘C. The past mean annual air temperature modelled experimentally for two sites in the Czech Republic hosting relict cryoturbation structures was between $-\mathrm{7.0}\pm \mathrm{1.9}$ and $-\mathrm{3.2}\pm \mathrm{1.5}$ ^∘C, which is well in line with earlier reconstructions utilizing various palaeo-archives. These initial results are promising and suggest that the model could become a useful tool for reconstructing Quaternary palaeo-environments across vast areas of mid-latitudes and low latitudes where relict periglacial assemblages frequently occur, but their full potential remains to be exploited.

1 Introduction

Many mid- and low-latitude regions of the world host a number of distinctive landforms and subsurface structures collectively termed relict periglacial features, such as various kinds of patterned ground, cryoturbations, frost wedges, solifluction structures, and blockfields, which developed during cold periods of the Quaternary due to intense freeze–thaw activity taking place under seasonal frost or permafrost conditions. Commonly, these features occur in places where other palaeo-indicators are rare, absent, or emerged at different times, which enhances their relevance as archives of past environmental conditions (Washburn, 1980; French, 2017; Ballantyne, 2018). So far, relict periglacial features have been used to reconstruct past climates almost exclusively on the basis of the climate thresholds of their active counterparts that are mostly situated in present-day high-latitude periglacial environments (Ballantyne and Harris, 1994). Such empirical interpretations are, however, problematic because suitable analogues for past periglacial environments are rare, particularly due to substantial latitudinal contrasts in solar insolation (Williams, 1975; French, 2017), and even if they can be found, active periglacial features that occur there may have developed under different climate conditions than those prevailing at the present time (Uxa et al., 2017; Ballantyne, 2018). Climate controls on the development of most periglacial features are thus poorly known, usually implying broad ranges of climate conditions (Washburn, 1980; Harris, 1982, 1994; Karte, 1983; Wayne, 1983; Ballantyne and Harris, 1994; Huijzer and Isarin, 1997; Ballantyne, 2018), which also relates to the fact that the features partly depend on other factors, such as ground physical properties, hydrology, topography, and ground-surface cover (Ballantyne, 2018). Consequently, the inferred palaeo-climates have frequently been thought to be of limited validity, and indeed most periglacial features have been widely accepted only as indicators of seasonal frost or permafrost and ground-ice presence (Ballantyne and Harris, 1994; Ballantyne, 2018). This situation can be largely attributed to prevailing interest in mapping the distribution patterns of periglacial features and their associations with mean annual air temperature (MAAT), which has been characteristic for traditional palaeo-periglacial geomorphology, while other details on their surface and subsurface dimensions have been widely overlooked. Greater emphasis on the surface and subsurface attributes of periglacial features, which are closely related to their formation and responsible processes, could, however, advance the discipline far beyond its current frontiers (see Barsch, 1993; French and Thorn, 2006).

Periglacial features form through various thermally induced and gravity-induced processes that mostly operate within a layer of seasonal freezing and thawing, the base of which is commonly sharply defined and confines the subsurface dimensions of the features (Williams, 1961). This zone is thus usually discernible in vertical cross sections because intense ice segregation and mass displacements associated with the formation of periglacial features alter the freeze–thaw layer so that its composition and properties differ from those of the underlying ground. This contrast may be preserved long after the periglacial features have ceased to be active, and, as such, it can indicate the thickness of the palaeo-freeze–thaw layer (French, 2017). Since the freeze–thaw depth is closely coupled with ground and air temperature conditions (e.g. Frauenfeld et al., 2004; Åkerman and Johansson, 2008; Wu and Zhang, 2010), it retains a valuable palaeo-climate record that can be approximated based on modern air temperature–freeze–thaw depth relations (Williams, 1975) or retrieved through an inverse solution of the equations calculating the freeze–thaw depth (Maarleveld, 1976; French, 2008). Obviously, this idea is not new, but despite its ingenuity, simplicity, and general acceptance in benchmark periglacial literature (e.g. Washburn, 1979; Ballantyne and Harris, 1994; French, 2017; Ballantyne, 2018), it has never been developed into a viable tool for deriving past thermal regimes that has a sound mathematical basis, is replicable, and lacks subjectivity (see Williams, 1975; Maarleveld, 1976) because computational methods have been durably underused by periglacial geomorphologists interested in reconstructions of Quaternary palaeo-environments.

This study presents and evaluates a simple modelling scheme called PERICLIMv1.0 (PERIglacial CLIMate) that is designed to infer palaeo-air temperature characteristics associated with relict periglacial features indicative of the palaeo-active-layer thickness, and the study discusses its uncertainties and applicability with respect to other palaeo-proxy records and/or model products. It specifically targets palaeo-active-layer phenomena because their palaeo-environmental significance and preservation potential are substantially higher than for seasonal frost features. Also, it intends to stimulate the application of modelling tools and foster the development of new quantitative methods in palaeo-environmental reconstructions utilizing relict periglacial features in order to improve their reputation as palaeo-proxy indicators.

2 Model description

The PERICLIMv1.0 builds on an inverse solution of the Stefan (1891) equation, which was originally developed to determine the thickness of sea ice, but it also describes the thaw propagation in ice-bearing grounds well, and lately it has become probably the most commonly used analytical tool for estimating the thickness of the active layer over permafrost (e.g. Klene et al., 2001; Shiklomanov and Nelson, 2002; Hrbáček and Uxa, 2020). It assumes that the thawed-zone temperature, which is controlled by the ground-surface temperature at the surface boundary, decreases linearly towards the bottom frozen zone that is constantly at 0 ^∘C and that latent heat is the only energy sink associated with its thawing, meaning that heat conduction below the thaw front is not accounted for (Kurylyk, 2015). As such, the Stefan equation tends to deviate inversely proportional to the moisture content in the active layer and the active-layer temperature at the onset of thawing (Romanovsky and Osterkamp, 1997; Kurylyk and Hayashi, 2016), but its accuracy is still reasonable (e.g. Klene et al., 2001; Shiklomanov and Nelson, 2002; Hrbáček and Uxa, 2020). Also, its simplicity and low requirements for input data compared to more complex analytical or numerical models (e.g. Kudryavtsev et al., 1977; Jafarov et al., 2012; Westermann et al., 2016) are highly advantageous for palaeo-applications as fewer assumptions have to be made. Here, it is solved for uniform, non-layered ground while ignoring any of its thaw-related mechanical responses.

The PERICLIMv1.0 deduces the thawing season temperature conditions in the above way, which it further converts into annual and freezing season air temperature attributes based on the assumed annual air temperature range as detailed below. Note that its code is implemented and disseminated as an R package.

Table 1 List of model input and output variables, their symbols, value ranges, and units.

2.1 Model input variables

The model requires inputs on the palaeo-active-layer thickness ξ [m], the volumetric ground moisture content ϕ [−], the dry ground bulk density ρ [kg m⁻³], the ground quartz content q [−], the ground grain-size class $f/c$ [−], the ground-surface thawing n factor n_t [−], and the annual air temperature range A_a [^∘C] (Table 1), which corresponds to the difference between the mean air temperature of the warmest and coldest month. The ground physical properties are assumed to characterize the entire modelling domain (∼active layer) and are treated as constant over time.

2.2 Ground-surface and air thawing index

The Stefan equation for calculating the active-layer thickness in a homogeneous substratum with constant physical properties has the following form (Lunardini, 1981):

$$\begin{array}{}\text{(1)}& \mathit{\xi }=\sqrt{{\displaystyle \frac{\mathrm{2}{k}_{\mathrm{t}}{I}_{\text{ts}}}{L\mathit{\varphi }{\mathit{\rho }}_{\mathrm{w}}}}},\end{array}$$

where k_t [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] is the thermal conductivity of the thawed ground calculated here as a function of its dry bulk density and volumetric moisture content using the Johansen (1977) thermal conductivity model (Appendix A), I_ts [^∘C d] is the ground-surface thawing index defined as a sum of positive daily ground-surface temperatures in the thawing season (Fig. 1), L [334 000 J kg⁻¹] is the specific latent heat of fusion of water, and ρ_w [1000 kg m⁻³] is the density of water. Note that I_ts must be multiplied by the scaling factor of 86 400 s d⁻¹ in Eq. (1) to obtain the active-layer thickness in metres. Also, the product of ϕ and ρ_w can be alternatively substituted by that of the gravimetric ground moisture content and dry ground bulk density as their results are identical.

Figure 1 An idealized course of air temperature during the year described by a sine function with a mean of −4 ^∘C and a range of 20 ^∘C (black solid line) superimposed on the annual air temperature curve with daily variations (grey solid line). The air thawing index is shown in red, while the blue areas depict partial air freezing indices for the preceding (left) and subsequent (right) freezing season, respectively. See the text and Table 1 for abbreviations.

The ground-surface thawing index required to reach a given active-layer thickness can be obtained if Eq. (1) is rearranged as

$$\begin{array}{}\text{(2)}& I_{\text{ts}}={\displaystyle \frac{{\mathit{\xi }}^{\mathrm{2}}L\mathit{\varphi }{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{2}{k}_{\mathrm{t}}}}.\end{array}$$

The ground-surface thawing index can then be converted into the air thawing index I_ta [^∘C d] through the ground-surface thawing n factor (Lunardini, 1978), which is a simple empirical transfer function that has been widely used to parametrize the thawing season air–ground temperature relations across permafrost landscapes (e.g. Klene et al., 2001; Gisnås et al., 2017):

$$\begin{array}{}\text{(3)}& I_{\text{ta}}={\displaystyle \frac{I_{\text{ts}}}{n_{\mathrm{t}}}}.\end{array}$$

Note that the effect of snow cover does not have to be accounted for because the Stefan equation (Eq. 2) considers solely the thawing season temperatures responsible for the active-layer thawing. Likewise, the resulting air thawing index is later used only to calculate other air temperature characteristics, and these are not affected by snow in any way. Such a scheme is particularly advantageous because it keeps the number of inputs small. Yet, it should be borne in mind that it is especially suitable for locations without long-lasting snow cover, which might disrupt the coupling between the air and ground temperatures during prolonged snow melting periods (e.g. Gisnås et al., 2016).

2.3 Air temperature characteristics

The temporal evolution of air temperature over the year T_a(t) [^∘C] can be described by a sine wave (Fig. 1) such as

$$\begin{array}{}\text{(4)}& T_{\mathrm{a}}\left(t\right)=\text{MAAT}+{\displaystyle \frac{A_{\mathrm{a}}}{\mathrm{2}}}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{P}}\right),\end{array}$$

where t [d] is the time and P [365 d] is the period of air temperature oscillations. Note that the annual air temperature range must be halved in Eq. (4) and the subsequent equations in order to characterize the annual temperature variations around MAAT (Fig. 1).

The air thawing index represents the positive area under the annual air temperature curve (Fig. 1) and can be calculated by integrating Eq. (4) over the thawing season as follows:

$$\begin{array}{}\text{(5)}& I_{\text{ta}}=\underset{t_{\mathrm{1}}}{\overset{t_{\mathrm{2}}}{\int }}{T}_{\mathrm{a}}\left(t\right)\mathrm{d}t,\end{array}$$

with

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{t}_{\mathrm{1}}=\arcsin \left(-{\displaystyle \frac{\text{MAAT}}{\frac{A_{\mathrm{a}}}{\mathrm{2}}}}\right){\displaystyle \frac{P}{\mathrm{2}\mathit{\pi }}},\text{(7)}& {\displaystyle }& {\displaystyle }{t}_{\mathrm{2}}=\left[\mathit{\pi }-\arcsin \left(-{\displaystyle \frac{\text{MAAT}}{\frac{A_{\mathrm{a}}}{\mathrm{2}}}}\right)\right]{\displaystyle \frac{P}{\mathrm{2}\mathit{\pi }}},\end{array}$$

where t₁ is the time when the air temperature curve crosses the 0 ^∘C level from below (∼thawing season begins), while t₂ is the time when it crosses this level from above (∼ thawing season ends) (e.g. Nelson and Outcalt, 1987). Unfortunately, Eq. (5) has no analytical solution for MAAT. However, it can be derived from a nomogram (Fig. 2) or, as here, it can be calculated numerically using the bisection root-finding method searching for MAAT such that $-\frac{A_{\mathrm{a}}}{\mathrm{2}}<\text{MAAT}\le \mathrm{0}$. This condition ensures that both positive and negative air temperatures have occurred during the year, which is an essential prerequisite for the active layer to form. Admittedly, it is simplistic because air–ground temperatures are modulated by surface and subsurface offsets so that the permafrost–seasonal frost boundary usually occurs at slightly negative MAAT (Smith and Riseborough, 2002). Consequently, there might be a risk that the model is incorrectly applied to seasonal frost conditions. However, this can be easily prevented if periglacial features that have indisputably developed in the presence of permafrost are examined.

Figure 2 A nomogram showing relations between the air thawing and freezing index, mean annual air temperature (solid diagonal lines), annual air temperature range (dashed curved lines), and mean air temperature of the warmest month (dotted curved lines). Note that the value of the mean annual air temperature (read diagonally) and the air freezing index (read horizontally) is obtained at the intersection of the air thawing index (read vertically) and the annual air temperature range or the mean air temperature of the warmest month (both read curvilinearly).

Once MAAT is known, the air freezing index I_fa [^∘C d] can be simply computed as

$$\begin{array}{}\text{(8)}& I_{\text{fa}}=\text{MAAT}P-I_{\text{ta}}.\end{array}$$

Furthermore, the mean air temperature of the warmest MATWM [^∘C] and coldest MATCM [^∘C] month is calculated as

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }\text{MATWM}=\text{MAAT}+{\displaystyle \frac{A_{\mathrm{a}}}{\mathrm{2}}},\text{(10)}& {\displaystyle }& {\displaystyle }\text{MATCM}=\text{MAAT}-{\displaystyle \frac{A_{\mathrm{a}}}{\mathrm{2}}}.\end{array}$$

The mean air temperature of the thawing MATTS [^∘C] and freezing MATFS [^∘C] season is defined as

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }\text{MATTS}={\displaystyle \frac{I_{\text{ta}}}{L_{\mathrm{t}}}},\text{(12)}& {\displaystyle }& {\displaystyle }\text{MATFS}={\displaystyle \frac{I_{\text{fa}}}{L_{\mathrm{f}}}},\end{array}$$

where L_t [d] and L_f [d] are the duration of the thawing and freezing season, respectively, which is expressed from Eqs. (6) and (7) as

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{t}}=t_{\mathrm{2}}-t_{\mathrm{1}}=\left[\mathit{\pi }-\mathrm{2}\arcsin \left(-{\displaystyle \frac{\text{MAAT}}{\frac{A_{\mathrm{a}}}{\mathrm{2}}}}\right)\right]{\displaystyle \frac{P}{\mathrm{2}\mathit{\pi }}},\text{(14)}& {\displaystyle }& {\displaystyle }{L}_{\mathrm{f}}=P-L_{\mathrm{t}}.\end{array}$$

Unconventionally, the model can also be solved with the range of annual air temperature oscillations defined by MATWM (see Williams, 1975) if its difference from MAAT (i.e. MATWM−MAAT) is substituted for $\frac{A_{\mathrm{a}}}{\mathrm{2}}$ in Eq. (4) and elsewhere (Appendix B). Unsurprisingly but importantly, solutions based on alternate model input variables, as suggested above, produce identical outputs, allowing model adaptations to specific situations and available data.

3 Model validation

The performance of the model was tested using the same simulation schemes as for palaeo-applications detailed in the next section (see Sect. 4), but on the basis of data from modern permafrost environments of the James Ross Island and the Alaskan Arctic (Appendix C). Comparisons of the modelled and observed data showed relatively good agreement and clustering along the lines of equality for most air temperature characteristics (Fig. 3), which suggests that the model might also work well over a wider range of climates. Generally, however, the model underestimated the results, with those from James Ross Island being somewhat more accurate and less scattered than those from Alaska (Fig. 3).

Figure 3 Means and standard deviations of observed and modelled air temperature characteristics at the James Ross Island and Alaskan Arctic validation sites. Acronyms ME and MAE under the plot labels correspond to the site-weighted mean error and the site-weighted mean absolute error, respectively, while those at the bottom right of the plots indicate the names of the validation sites: Abernethy Flats (AF), Berry Hill slopes (BHS), Johann Gregor Mendel (JGM), Johnson Mesa (JM), Deadhorse (DH), Franklin Bluffs (FB), and West Dock (WD).

MAAT exhibited an average mean error and an average mean absolute error of −1.3 and 1.3 ^∘C, respectively (Fig. 3). Mean absolute error was ≤1 and ≤2 ^∘C at 57 % and 86 % of the validation sites, respectively, and its maximum was 2.4 ^∘C.

MATWM and MATCM showed slightly lower bias as their average mean errors attained −1.0 and −0.9 ^∘C, respectively, and average mean absolute errors achieved 1.1 ^∘C for both characteristics (Fig. 3). Their mean absolute deviations were ≤1 ^∘C at 43 % of the validation sites and ≤2 ^∘C at 100 % of them. Mean absolute MATWM and MATCM errors reached a maximum of 1.9 and 1.8 ^∘C, respectively.

I_ta and I_fa were modelled with average mean errors of −60 and −403 ^∘C d, respectively, which corresponds to −16 % and −12 % of the observed average mean values, and their average mean absolute errors reached 71 and 403 ^∘C d, respectively (Fig. 3). I_ta was biased by ≤40 ^∘C d at 29 % of the validation sites and by ≤80 ^∘C d at 43 % of them. I_fa, which reaches 1 order of magnitude larger values, deviated by ≤400 ^∘C d at 57 % of the validation sites and by ≤800 ^∘C d at 100 % of them. Maximum mean absolute departures of I_ta and I_fa reached 106 and 789 ^∘C d, respectively.

MATTS exhibited an average mean error and an average mean absolute error of −0.5 and 0.5 ^∘C, respectively, while for MATFS the mean errors averaged −0.6 and 1.1 ^∘C (Fig. 3). The agreement between the modelled and observed MATTS was ≤1 ^∘C at 100 % of the validation sites, with a maximum mean absolute error of 1.0 ^∘C. MATFS deviated by ≤1 ^∘C at 71 % of the validation sites, and at worst it was 2.5 ^∘C.

Since L_t and L_f inherently counteract, their characteristics mirror each other if the values are not rounded. L_t tended to be underestimated by an average of 21 d, while L_f was overestimated by an average of 20 d (Fig. 3), which comprised −20 % and 8 % of the observed average mean L_t and L_f, respectively, and their average mean absolute errors achieved 27 and 26 d. Mean underestimation or overestimation was ≤20 d at 43 % of the validation sites and ≤40 d at 86 %. Maximum mean deviation was up to 49 d.

4 Model application

As a feasibility study, the model was utilized experimentally for the derivation of palaeo-air temperature conditions on the basis of relict cryoturbation structures. Cryoturbations (∼periglacial involutions) are characterized by folded and/or dislocated strata of unconsolidated sediments caused by recurrent freeze–thaw-induced processes operating within the active layer over permafrost, which limits the vertical extent of the cryoturbations from below and also acts as an impervious boundary that provokes well-saturated conditions. As such, cryoturbations are thought to indicate the thickness of the active layer and the presence of permafrost at the time of their development. Also, MAAT thresholds of $<-\mathrm{8}$ and $<-\mathrm{4}$ ^∘C have been suggested for their formation within coarse- and fine-grained substrates, respectively (Vandenberghe, 2013; French, 2017), the validity of which can be assessed with the model as well.

4.1 Study sites

We consider two study sites in the Czech Republic where vertical cross sections through cryoturbated horizons, portraying the palaeo-active layers (Fig. 4), were exposed and sampled for attributes that allowed us to define the most plausible ranges of the model input variables.

Figure 4 Cross sections through cryoturbation structures at the (a) Brno–Černovice and (b) Nebanice sites, with the white dashed lines indicating their vertical extent and thereby the thickness of the palaeo-active layer. Note that the sections above the suggested base of the palaeo-active layer show distinct involution structures caused by freeze–thaw-induced processes, while those below have largely retained their primary structures with horizontal to sub-horizontal bedding planes.

The Brno–Černovice site ($\mathrm{49}{}^{\circ }{\mathrm{10}}^{\prime }{\mathrm{43}}^{\prime \prime }$ N, $\mathrm{16}{}^{\circ }{\mathrm{38}}^{\prime }{\mathrm{56}}^{\prime \prime }$ E; 240 m a.s.l.) is an active sand–gravel pit situated about 3 km southeast of the centre of Brno, southern Moravia. The mine is embedded within sands and gravels of the Tuřany terrace of the Svitava River, in reality a highly flattened alluvial fan tentatively attributed to the Günz glaciation (∼Marine Isotope Stage 22–11), which is located about 1.5 km east of the present channelized riverbed and ∼40 m higher. The material tends to coarsen down to a depth of 6–13.5 m, under which the Neogene clayey sands and gravels emerge (Musil et al., 1996; Musil, 1997; Bubík et al., 2000). The cryoturbations mostly consist of sands to gravels that constitute ∼80 % and ∼14 % of the substrate, respectively, and are deformed by numerous injection tongues (Fig. 4). Currently, the cryoturbations rest ∼1.6–3.6 m under the ground surface, but we suppose they were just below it when they developed.

The Nebanice site ($\mathrm{50}{}^{\circ }{\mathrm{07}}^{\prime }{\mathrm{13}}^{\prime \prime }$ N, $\mathrm{12}{}^{\circ }{\mathrm{27}}^{\prime }{\mathrm{09}}^{\prime \prime }$ E; 430 m a.s.l.) is a currently inactive sand–gravel pit about 1.5 km west of the centre of Nebanice, western Bohemia. It is set in sands and gravels of the Nebanice terrace of the Ohře River, which is thought to have originated during the Mindel–Riss glaciations (∼ Marine Isotope Stage 10) and is now situated about 200 m north of the riverbed at a relative height of ∼10 m (Šantrůček et al., 1994; Balatka et al., 2019). The terrace sediments have a thickness of 1.7–3.8 m (Balatka et al., 2019) and are underlain by the Pliocene–lower Pleistocene clays, sands, and gravels (Špičáková et al., 2000). The cryoturbations are mostly composed of sands and gravels that comprise ∼67 % and ∼29 % of the material, respectively, and take the form of festoons and ball-and-pillow structures (Fig. 4) emerging in the uppermost part of the terrace ∼1.1–2.8 m below the modern ground surface as they were later overlaid by younger sediments and soils as well.

Cryoturbations, like other permafrost-related features, in the area of the Czech Republic have been tentatively attributed to the last glacial period (Czudek, 2005). We believe that the studied relict cryoturbations indicate past environmental conditions similar to those at the end of the Last Glacial Maximum (LGM) when analogous features formed in nearby regions (Kasse, 1993; Huijzer and Vandenberghe, 1998; Kasse et al., 2003; Bertran et al., 2014).

4.2 Model set-up

4.2.1 Palaeo-active-layer thickness

The palaeo-active-layer thickness was considered to be represented by the vertical extent of the cryoturbated horizons (Fig. 4), which was determined in horizontal steps of 0.2 m along the entire length of each cross section (4.4 and 8.6 m), giving 22 and 43 measurements that averaged 1.58±0.28 and 1.40±0.13 m at the Brno–Černovice and Nebanice sites, respectively (Table 2).

Table 2 Variables used to model palaeo-air temperature characteristics related to the cryoturbations at the Brno–Černovice and Nebanice sites. Note that the superscripts norm, beta, and unif indicate normal, beta, and uniform distributions describing the respective variables.

See Table 1 for abbreviations.

4.2.2 Ground physical properties

Ground physical properties were set using intact samples collected at four representative positions along four vertical profiles within each cross section (16 samples per cross section) into 100 cm³ stainless-steel cylinders, which were then weighted in wet and dry states as well as sieved to assess their volumetric moisture content, dry bulk density, and texture. The current volumetric moisture content averaged 8.8 % and 11.4 % at the Brno–Černovice and Nebanice sites, respectively, and these values were assumed to be the lower moisture thresholds. The upper moisture thresholds were supposed to be given by average ground porosity, which was calculated as a function of the dry bulk density (Eq. A4), and reached 39.4 % and 39.1 %, respectively (Table 2). Because cryoturbations require well-saturated conditions for at least part of the thawing season for their development (Vandenberghe, 2013; French, 2017), we further skewed the moisture content values leftwards using a beta distribution as follows:

$$\begin{array}{}\text{(15)}& \mathit{\varphi }\left(x\right)={\mathit{\varphi }}_{\mathrm{0}}+f(x;\mathit{\alpha },\mathit{\beta })\left(n-{\mathit{\varphi }}_{\mathrm{0}}\right),\end{array}$$

where ϕ₀ [−] is the current average volumetric ground moisture content, $f(x;\mathit{\alpha },\mathit{\beta })$ expresses the probability density function of the beta distribution for $\mathrm{0}\le x\le \mathrm{1}$ with shape parameters tentatively assumed to be α=5 and β=2 that yield a mode of 0.8, and n [−] is the average ground porosity. The resulting values are thus bounded by the lower and upper moisture thresholds (Table 2) but show left-skewed distributions peaking at ∼33.3 % and ∼33.6 % at the Brno–Černovice and Nebanice sites, respectively, which corresponds to a degree of saturation of ∼84.5 % and ∼85.8 % (Fig. 5). Since the vast majority of moisture in sands and gravels tends to undergo phase changes (Andersland and Ladanyi, 2004) and its unfreezing portion is supposed to influence the calculations negligibly if its ratio to the total moisture content is at levels up to several tens of percent (Uxa, 2017), the moisture contents were not further adjusted for unfrozen moisture. Slightly lower moisture content of the eastern site in Brno–Černovice compared to Nebanice in the west (Fig. 5) is likely to be reasonable because it can be understood as including the continentality and elevation effects as well as the thicker palaeo-active layer (Table 2) in which the same amount of water has a lower volume fraction.

Figure 5 Probability distributions of volumetric ground moisture contents assumed for the Brno–Černovice (BC) and Nebanice (NB) sites.

Other ground physical properties were supposed to be unchanged since the cryoturbations developed. The dry ground bulk density was 1635±120 and 1645±116 kg m⁻³ at the Brno–Černovice and Nebanice sites, respectively. The ground quartz content was estimated at 30 %–56 % and 30 %–57 %, respectively, on the basis of the proportions of clay–silt and sand fractions (see Appendix A). Given the low proportion of clay–silt fraction and the presence of gravel (see Sect. 4.1), the substrates were treated as coarse-grained (sensu Johansen, 1977) at both sites (Table 2).

4.2.3 Ground-surface thawing n factor

Since the coupling between the air and ground-surface temperatures is principally governed by ground-surface cover and its properties themselves (Westermann et al., 2015), the ground-surface thawing n factors were estimated on the basis of values published elsewhere for bare (excluding bedrock and debris) to sparsely vegetated (lichens, mosses, or grasses) surfaces, which are assumed to have existed at the study sites when the cryoturbations originated because treeless landscapes dominated in southern Moravia and western Bohemia at that time (e.g. Kuneš et al., 2008) and also because the cryoturbated horizons contain no or negligible amounts of organic remains. We took a total of 41 ground-surface thawing n-factor values reported from mid-latitude and mostly permafrost regions of central–eastern Norway (Juliussen and Humlum, 2007) and British Columbia–Yukon and Labrador, Canada (Lewkowicz et al., 2012; Way and Lewkowicz, 2018), that are believed to be representative for the study sites. Those locations also occur outside the Arctic Polar Circle and, as such, have comparable insolation budgets owing to the absence of polar-day periods, which might impose an undesirable growth of the ground-surface thawing n-factor values (e.g. Shur and Slavin-Borovskiy, 1993). Overall, the collected n factors averaged 1.03±0.12, which was assigned to both study sites (Table 2).

4.2.4 Annual air temperature range

As a starting point, we used the annual air temperature ranges based on monthly air temperatures measured in 1981–2010 at meteorological stations located about 4 and 7 km southeast and southwest of the Brno–Černovice and Nebanice sites, respectively, at elevations of 241 and 475 m a.s.l. (Czech Hydrometeorological Institute, 2020), which averaged 23.2±2.4 and 20.9±2.6 ^∘C, respectively (Table 2). Additionally, we also assumed stepwise perturbations of 2 ^∘C to these annual air temperature ranges up to 33.2±2.4 and 30.9±2.6 ^∘C, respectively, which is within the range of 28–36 ^∘C suggested by previous regional estimates for the end of the LGM (Huijzer and Vandenberghe, 1998) and at the same time maintains the contrasts between the study sites (Table 2).

4.2.5 Simulations

The model was solved by the Monte Carlo simulation with Latin hypercube sampling (LHS) (McKay et al., 1979) using the lhs R package (Carnell, 2020), which discretizes the probability density functions of the individual model input variables (Table 2) into N non-overlapping intervals of equal probability $\mathrm{1}/N$ and then randomly selects one value from each. Subsequently, these samples are randomly matched and used as uncorrelated inputs for N model runs (McKay et al., 1979). LHS is computationally more efficient than simple random sampling because its sampling scheme adapts to the probability density functions of the input variables, and thus it requires fewer simulations to achieve stable outputs if N is large enough.

We realized 1000 model runs for six scenarios of the annual air temperature ranges at each study site, of which 78.1 %–91.1 % produced physically feasible combinations of the input variables and provided the most plausible palaeo-air temperature characteristics that account for the uncertainty and natural temporal variability of the inputs.

4.3 Modelled palaeo-air temperature characteristics

Most modelled palaeo-air temperature characteristics changed to various extents depending on the scenarios of the annual air temperature ranges, the higher values of which generally caused colder and/or more continental conditions as well as larger scatter in the model outputs. This was especially true for MAAT and the freezing season characteristics, whereas the changes were comparatively milder for the thawing season attributes (Fig. 6).

Figure 6 Means and standard deviations of modelled palaeo-air temperature characteristics at the Brno–Černovice (BC) and Nebanice (NB) sites for six scenarios of the annual air temperature ranges.

MAAT was modelled at $-\mathrm{6.6}\pm \mathrm{2.7}$ to $-\mathrm{3.3}\pm \mathrm{2.0}$ ^∘C and $-\mathrm{7.0}\pm \mathrm{1.9}$ to $-\mathrm{3.2}\pm \mathrm{1.5}$ ^∘C at the Brno–Černovice and Nebanice sites, respectively (Fig. 6), which corresponds to its average decline of −16.0 to −12.7 ^∘C and −15.1 to −11.3 ^∘C, respectively, in comparison with the 1981–2010 period. The average contrasts between the warmest and coldest MAAT scenarios based on the current and assumed past annual air temperature ranges were 3.3 and 3.8 ^∘C at the Brno–Černovice and Nebanice sites, respectively, which means that MAAT changed by 33 % and 38 % of the change in the annual air temperature range.

The thawing seasons had MATTS of 5.4±1.2 to 6.5±1.6 ^∘C and 4.7±0.8 to 5.5±1.0 ^∘C at the Brno–Černovice and Nebanice sites, respectively, and culminated with MATWM of 8.4±1.9 to 10.1±2.6 ^∘C and 7.3±1.3 to 8.5±1.6 ^∘C (Fig. 6), which suggests that MATTS changed on average by 11 % and 8 % of the magnitudes of the annual air temperature range perturbations, while MATWM varied on average by 17 % and 12 %. I_ta showed even more consistent outputs as it reached 823±271 to 915±353 ^∘C d and 704±181 to 721±203 ^∘C d at the Brno–Černovice and Nebanice sites, respectively (Fig. 6). As such, it varied by as little as ∼1.1 % and ∼0.2 % per 1 ^∘C change in the annual air temperature range. Likewise, L_t was modelled at 135±20 to 149±20 d and 128±14 to 147±15 d (Fig. 6), which yields respective changes of ∼1.0 % and ∼1.5 %.

The freezing seasons exhibited the largest scatter and the highest variability among the scenarios, with MATFS of $-\mathrm{14.2}\pm \mathrm{1.8}$ to $-\mathrm{9.3}\pm \mathrm{1.5}$ ^∘C and $-\mathrm{13.7}\pm \mathrm{1.6}$ to $-\mathrm{8.5}\pm \mathrm{1.5}$ ^∘C at the Brno–Černovice and Nebanice sites, respectively, and MATCM as low as $-\mathrm{23.2}\pm \mathrm{3.2}$ to $-\mathrm{15.0}\pm \mathrm{2.6}$ ^∘C and $-\mathrm{22.5}\pm \mathrm{2.8}$ to $-\mathrm{13.7}\pm \mathrm{2.5}$ ^∘C (Fig. 6). It thus follows that their variations reached on average 49 % and 52 % of the annual air temperature range perturbations for MATFS and as much as 82 % and 88 % for MATCM. I_fa spanned $-\mathrm{3309}\pm \mathrm{667}$ to $-\mathrm{2041}\pm \mathrm{492}$ ^∘C d and $-\mathrm{3270}\pm \mathrm{523}$ to $-\mathrm{1873}\pm \mathrm{429}$ ^∘C d at the Brno–Černovice and Nebanice sites, respectively (Fig. 6), which corresponds to average variations of ∼6.2 % and ∼7.5 % per 1 ^∘C change in the annual air temperature range. On the other hand, L_f was modelled at 216±20 to 230±20 d and 218±15 to 237±14 d, respectively, and thus it varied on average by as little as ∼0.6 % and ∼0.9 % of the magnitudes of the annual air temperature range perturbations.

5 Discussion

5.1 Model performance and limitations

Generally, the model validation using the modern data (Sect. 3) showed relatively high accuracy for most air temperature characteristics if suitable input variables were available. It is remarkable, though, given that active-layer thickness commonly exhibits rather moderate to low correlations with annual or freezing season air and ground temperature attributes, but on the other hand, it mostly strongly couples with thawing season air and ground temperature indices (e.g. Frauenfeld et al., 2004; Åkerman and Johansson, 2008; Wu and Zhang, 2010). Since the model builds on active-layer thickness–thawing season temperature relations, which it further converts into annual and freezing season air temperature characteristics, this scheme gives rise to reasonable accuracy. It also needs to be stressed that the model was successful regardless of the fact that the ground physical properties used have certainly undergone at least slight changes since sampling and varied over the validation period. More accurate and less scattered outputs at the James Ross Island sites were largely due to a rather homogeneous distribution of ground physical properties within the active layer there (Hrbáček et al., 2017a; Hrbáček and Uxa, 2020). By contrast, the Alaskan profiles have a two-layer composition with peat over mineral soil (Zhang, 1993). Additionally, the model parameterizes the air temperature behaviour with a sine wave, which simplifies its actual evolution over the year and completely ignores sub-annual variations.

Overall, however, the model underestimated most air temperature characteristics (Fig. 3). The underestimation is attributed to intrinsic shortcomings of the Stefan equation (Eq. 2) that tends to deviate inversely proportional to the moisture content in the active layer and the active-layer temperature at the onset of thawing (Romanovsky and Osterkamp, 1997; Kurylyk and Hayashi, 2016). Also, it is associated with the Johansen (1977) thermal conductivity model (Appendix A), which tended to produce overly high thermal conductivity values at the James Ross Island sites. Surely, the Stefan equation (Eq. 2) might be improved by a number of correction factors, but these require additional inputs, such as frozen thermal conductivity, thawed and frozen volumetric heat capacity, or active-layer temperature at the start of thawing (Kurylyk and Hayashi, 2016). As such, the corrections are frequently difficult to implement even in many present-day situations and are definitely much less viable for palaeo-applications. Moreover, their inverse solution would not be straightforward and would probably demand iterative techniques. Similarly, the Johansen (1977) thermal conductivity model is advantageous in that it requires fewer inputs compared with other solutions while having comparable accuracy (e.g. Dong et al., 2015; He et al., 2017; Zhang et al., 2018). As such, it is also difficult to replace by other thermal conductivity schemes.

It should also be highlighted that for the model validation cases the active-layer thickness was reduced by the depth of the temperature sensor, which was used to determine the ground-surface thawing n factor (see Appendix C). However, if this treatment was not done, the model outputs improved, with MAAT, MATWM, MATCM, MATTS, and MATFS showing average errors ≤1 ^∘C, I_ta and I_fa deviating on average by −1 % and −8 %, respectively, and L_t and L_f by −16 % and 6 %, respectively, because this compensated for the intrinsic deviations of the Stefan equation (Eq. 2) described in the previous paragraph. Since the published ground-surface thawing n factors used for the palaeo-applications (Juliussen and Humlum, 2007; Lewkowicz et al., 2012; Way and Lewkowicz, 2018) built on various near-surface depths of ground temperature sensors, we did not adjust the active-layer thickness in these cases, and the modelled palaeo-air temperature characteristics could thus have slightly increased accuracy as well.

5.2 Comparison to previous palaeo-air temperature reconstructions

The palaeo-MAAT modelled for two sites in the Czech Republic (Sect. 4.3) was between $-\mathrm{7.0}\pm \mathrm{1.9}$ and $-\mathrm{3.2}\pm \mathrm{1.5}$ ^∘C and its corresponding reduction between −16.0 and −11.3 ^∘C in comparison with the 1981–2010 period, which is relatively consistent with earlier reconstructions utilizing various relict periglacial features in the central European lowlands that suggested MAAT depressions mostly between −17 and −12 ^∘C (Poser, 1948; Büdel, 1953; Kaiser, 1960; Frenzel, 1967; Goździk, 1973; Huijzer and Vandenberghe, 1998; Marks et al., 2016). By contrast, it disagrees with slightly milder MAAT reductions of at least −7 and −13 to −6 ^∘C derived from groundwater data (Corcho Alvarado et al., 2011) and borehole temperature logs (Šafanda and Rajver, 2001), but these may not necessarily correspond to the lowest temperatures because groundwater cycling has been slowed or interrupted by permafrost, while ground temperature history may have been partly masked by latent heat effects. Similarly, the modelled MAAT differs rather highly from simulations of three global climate models (GCMs), namely the Community Climate System Model 4, the Model for Interdisciplinary Research on Climate Earth System Model, and the Max Planck Institute Earth System Model Paleo, capturing the LGM at ∼22 ka, which were originally released by the Coupled Model Intercomparison Project 5 in coarse resolutions but later downscaled to ${\mathrm{2}}^{\prime }{\mathrm{30}}^{\prime \prime }$ and made available via the WorldClim 1.4 dataset at https://www.worldclim.org (last access: 8 August 2016) (Hijmans et al., 2005). The GCMs suggest that MAAT was between −4.5 and −0.4 ^∘C at the study sites, which corresponds to its reduction of between −12.6 and −9 ^∘C in comparison with the 1981–2010 period. Such relatively high temperatures, however, contrast with many proxy records and are thus suspected to be unrepresentative of the coldest LGM conditions when continuous permafrost presumably occurred there (Vandenberghe et al., 2014; Lindgren et al., 2016).

The modelled MATWM of 7.3±1.3 to 10.1±2.6 ^∘C is well in the range of proxy-based MATWM of 5 to 13 ^∘C that has been reconstructed for the central European lowlands (Huijzer and Vandenberghe, 1998; Marks et al., 2016). By contrast, it deviates greatly from the GCM outputs, being as high as 11.1 to 18 ^∘C. The modelled MATCM ranged widely from $-\mathrm{23.2}\pm \mathrm{3.2}$ to $-\mathrm{13.7}\pm \mathrm{2.5}$ ^∘C, which, however, also pertains to other proxy records that show somewhat lower values of −27 to −16.5 ^∘C (Huijzer and Vandenberghe, 1998; Marks et al., 2016). The GCM-based MATCM exhibits less variability of −21.6 to −17.7 ^∘C, but its average is generally consistent with that modelled by this study. Unfortunately, other modelled palaeo-air temperature characteristics cannot be validated directly because no reliable proxy records or GCM outputs are available for them. However, as they are closely related to MAAT, MATWM, and MATCM, we believe that those attributes have similar plausibility.

Lastly, the modelled MAAT between $-\mathrm{7.0}\pm \mathrm{1.9}$ and $-\mathrm{3.2}\pm \mathrm{1.5}$ ^∘C is relatively consistent with the MAAT threshold of $<-\mathrm{4}$ ^∘C commonly suggested for the formation of cryoturbation structures within fine-grained substrates (Vandenberghe, 2013). At the study sites, however, the cryoturbations consist of coarse-grained materials, for which a MAAT threshold as low as $<-\mathrm{8}$ ^∘C has been proposed (Vandenberghe, 2013). This study thus raises questions about the validity of the previously suggested MAAT thresholds for cryoturbation structures (see Vandenberghe, 2013; French, 2017) and calls for their thorough revision.

5.3 Sensitivity analysis

Global sensitivity analysis using multiple regression suggested that the palaeo-active-layer thickness and annual air temperature range had a major impact on the modelled palaeo-air temperature characteristics at the Brno–Černovice and Nebanice sites (Fig. 7). The palaeo-active-layer thickness importantly showed the highest values of the standardized regression coefficients (SRCs), especially for the annual and thawing season air temperature attributes. Similarly, the annual air temperature range also highly influenced MAAT, and, in particular, it had the utmost control over the freezing season air temperatures as its SRCs even tended to outweigh those for the palaeo-active-layer thickness at that time of the year. It thus follows that the freezing season characteristics may have limited accuracy if the annual air temperature range is uncertain, which indeed translates into their higher variance (Fig. 6). On the other hand, the annual air temperature range slightly affected the thawing season air temperature attributes (Fig. 7), and these are thus assumed to be the most plausible. Ground-surface and subsurface input variables, such as volumetric ground moisture content, ground dry bulk density, and ground-surface thawing n factor, had considerably lower, albeit stable, influences on most modelled palaeo-air temperature characteristics. Ground quartz content was the weakest of the input variables as it was responsible for only minor variability in the model outputs (Fig. 7).

Figure 7 Sensitivity of the modelled palaeo-air temperature characteristics at the Brno–Černovice (BC) and Nebanice (NB) sites to individual model input variables expressed by standardized regression coefficients. See the text and Table 1 for abbreviations.

5.4 Model applicability to periglacial features

Besides cryoturbations, we assume that the model could be utilized to derive palaeo-air temperature characteristics on the basis of any relict periglacial features that can indicate former active-layer thickness, such as some kinds of patterned ground, some solifluction structures, frost-wedge tops, autochthonous blockfields, mountaintop detritus, active-layer detachment slides, up-frozen clasts, indurated horizons, and frost weathering microstructures (Ballantyne and Harris, 1994; Matsuoka, 2011; Ballantyne, 2018). Also, it could be easily adapted for seasonal frost features, although the estimation of snow conditions would be complicated.

Since most periglacial features develop on at least decadal to centennial timescales (e.g. Karte, 1983; Matsuoka, 2001; Ballantyne, 2018) inherently involving natural variations in climate and active-layer thickness, we hypothesize that their vertical extent frequently also includes the bottom transient layer in which the boundary between the active layer and permafrost fluctuates when periglacial features form (see Shur et al., 2005). This implies that the palaeo-active-layer thickness may appear as a dispersed rather than a sharp boundary in many vertical cross sections. Special attention must thus be paid to avoid any ambiguity in both the identification of relict periglacial features and the determination of the palaeo-active layer that may be severely degraded. Indeed, some periglacial features may be produced by seasonal frost alone, while some lookalike features may even have a non-periglacial origin (Ballantyne and Harris, 1994; Ballantyne, 2018). Nonetheless, even if identified correctly, problems may arise in places where the ground-surface level has changed over time due to sedimentation or erosion. High uncertainty can probably be expected for periglacial features composed of pebbly to bouldery materials, such as blockfields or mountaintop detritus, in which the vertical variability of ground physical properties is extremely high (Ballantyne, 1998) and difficult to describe by single-value variables. On slopes, these materials may also provoke non-conductive heat transfer processes, which give rise to high-magnitude variations in ground temperatures over short distances that cannot be addressed by simple heat conduction models (Wicky and Hauck, 2017).

Given the above considerations, modelling should ideally utilize coexisting periglacial features for more robust palaeo-air temperature estimates. Admittedly, it can capture only a snapshot of the temperature history, but this is no different from numerous other palaeo-indicators, such as various glacial deposits. Moreover, if relict periglacial assemblages of different ages occur in a given region, they may eventually provide a more complete record of former temperature conditions. Undoubtedly, dating of periglacial features is still challenging because they can have a highly complex formation history, which partly devaluates their palaeo-environmental importance. Nonetheless, this shortcoming is increasingly being suppressed by improved dating methods that bring more reliable periglacial chronologies (e.g. Andrieux et al., 2018; Nyland et al., 2020; Engel et al., 2021).

6 Conclusions

The PERICLIMv1.0 is a novel easy-to-use model that derives palaeo-air temperature characteristics related to the palaeo-active-layer thickness, which can be recognized using many relict periglacial features found in past permafrost regions. The model evaluation against modern temperature records demonstrated that it reproduces air temperature characteristics, such as MAAT, MATWM, MATCM, MATTS, and MATFS, with average errors ≤1.3 ^∘C. I_ta and I_fa deviate on average by −16 % and −12 %, respectively, while L_t and L_f tend to be on average underestimated and overestimated by −20 % and 8 %, respectively. The palaeo-MAAT modelled for two sites in the Czech Republic hosting relict cryoturbation structures was between $-\mathrm{7.0}\pm \mathrm{1.9}$ and $-\mathrm{3.2}\pm \mathrm{1.5}$ ^∘C, and its corresponding reduction was between −16.0 and −11.3 ^∘C in comparison with the 1981–2010 period, which is relatively well in line with earlier reconstructions utilizing various palaeo-archives.

These initial results are promising and suggest that the model could become a useful tool for reconstructing Quaternary palaeo-environments across vast areas of mid-latitudes and low latitudes where relict periglacial assemblages frequently occur, but their full potential remains to be exploited. It is the very first viable solution that seeks to interpret relict periglacial features quantitatively and in a replicable and subjectivity-suppressed manner. As such, it can provide much more plausible periglacial-based palaeo-air temperature reconstructions than before. Hopefully, it will be a springboard for follow-up developments of more sophisticated modelling tools that will further increase the exploitability and reliability of relict periglacial features as indicators of palaeo-climates.

Appendix A: Thermal conductivity of the thawed ground

Thermal conductivity of the thawed ground is calculated as a function of its dry bulk density and volumetric moisture content using the Johansen (1977) thermal conductivity model. This empirical transfer scheme requires fewer parameterizations than its derivatives, but it still provides reasonable and consistent estimates of thermal conductivity for a wide range of substrates having a degree of saturation of >5 %–10 % to 100 % (e.g. Dong et al., 2015; He et al., 2017; Zhang et al., 2018). Basically, it interpolates between the dry k_dry [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] and saturated k_sat [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] thermal conductivity of a material as follows:

$$\begin{array}{}\text{(A1)}& k_{\mathrm{t}}=k_{\text{dry}}+(k_{\text{sat}}-k_{\text{dry}})K_{\mathrm{e}},\end{array}$$

where K_e [−] is the dimensionless Kersten number used to normalize the contrast between k_sat and k_dry based on the degree of saturation.

The dry ground thermal conductivity is defined by the following semi-empirical relationship (Johansen, 1977):

$$\begin{array}{}\text{(A2)}& k_{\text{dry}}={\displaystyle \frac{\mathrm{0.135}\mathit{\rho }+\mathrm{64.7}}{\mathrm{2700}-\mathrm{0.947}\mathit{\rho }}},\end{array}$$

where the constant 2700 kg m⁻³ represents the typical density of solid ground particles that is used throughout the thermal conductivity scheme to ensure consistency of the calculations.

The saturated ground thermal conductivity is calculated by a weighted geometric mean based on the thermal conductivities of individual ground constituents and their respective volume fractions (Johansen, 1977):

$$\begin{array}{}\text{(A3)}& k_{\text{sat}}=k_{\mathrm{s}}^{\mathrm{1}-n}{k}_{\mathrm{w}}^n,\end{array}$$

where k_s [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] and k_w [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] are the thermal conductivity of solid ground particles and water, respectively, and n [−] is the ground porosity, which is expressed as a function of the dry bulk density and the typical density of solids:

$$\begin{array}{}\text{(A4)}& n=\mathrm{1}-{\displaystyle \frac{\mathit{\rho }}{\mathrm{2700}}}.\end{array}$$

The thermal conductivity of water is set at 0.57 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, while that of solids is computed as

$$\begin{array}{}\text{(A5)}& k_{\mathrm{s}}=k_{\mathrm{q}}^q{k}_{\mathrm{o}}^{\mathrm{1}-q},\end{array}$$

where k_q [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] and k_o [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] are the thermal conductivity of quartz and other minerals, respectively, and q [−] is the quartz fraction of the total content of solids. Quartz is assigned the thermal conductivity of 7.7 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, while for other minerals it is as follows (Johansen, 1977):

$$\begin{array}{}\text{(A6)}& k_{\mathrm{o}}=\left\{\begin{array}{l}\mathrm{3},\text{for}q<\mathrm{20}\mathit{\%}\wedge \text{coarse-grained}\\ \mathrm{2},\text{otherwise}\end{array}\right..\end{array}$$

Note that materials having more than 5 % clay should be considered fine-grained, while others should be treated as coarse-grained (sensu Johansen, 1977).

Since the quartz content is usually unknown, it can be estimated as a weighted average of its empirically obtained percentages for clay, silt, and sand fraction, having 0 %, 15 %, and 45 %, respectively, which is thought to have an uncertainty <30 % (Johansen, 1977). Alternatively, it can also be drawn from a nomogram using the ground texture (Fig. A1).

Figure A1 A nomogram for estimating the quartz content using the ground texture based on Johansen (1977). Note that the quartz content (black solid diagonal lines) is obtained at the intersection of the contents of clay, silt, and sand (grey dashed three-way lines) that are read in the direction of their respective tick marks and labels.

The Kersten number for thawed fine- and coarse-grained ground with a degree of saturation S [−] larger than 10 % and 5 %, respectively, is expressed as (Johansen, 1977)

$$\begin{array}{}\text{(A7)}& K_{\mathrm{e}}=\left\{\begin{array}{l}\log S+\mathrm{1},\text{for}S>\mathrm{10}\mathit{\%}\wedge \text{fine-grained}\\ \mathrm{0.7}\log S+\mathrm{1},\text{for}S>\mathrm{5}\mathit{\%}\wedge \text{coarse-grained}\end{array}\right.,\end{array}$$

where

$$\begin{array}{}\text{(A8)}& S={\displaystyle \frac{\mathit{\varphi }}{n}}.\end{array}$$

Clearly, the values of the calculated ground thermal conductivity increase exponentially and logarithmically with rising input values of the dry bulk density and volumetric moisture content, respectively, and thus the conductivity tends to change more sharply if the inputs are higher and lower, respectively (Fig. A2).

Figure A2 Thermal conductivity of the thawed ground as a function of its dry bulk density and volumetric moisture content calculated using the Johansen (1977) thermal conductivity model for thawed fine substrates with a degree of saturation of >10 % to 100 % and a quartz content of 40 %. Note that the white areas are outside the saturation limits.

Appendix B: Solution using MATWM to define the range of annual air temperature oscillations

Unconventionally, the range of annual air temperature oscillations can also be defined using MATWM (see Williams, 1975) if its difference from MAAT (i.e. MATWM−MAAT) is substituted for $\frac{A_{\mathrm{a}}}{\mathrm{2}}$ in Eq. (4) and elsewhere. Likewise, MAAT is then calculated numerically using the bisection root-finding algorithm, but it is searched for $-\mathrm{\infty }<\text{MAAT}\le \mathrm{0}$ because the actual value of the annual air temperature range is to be determined by the calculation itself. This solution can be advantageously combined with other palaeo-indicators that allow for the estimation of MATWM. However, it should be employed cautiously because it is highly sensitive to variations of I_ta and MATWM itself (Fig. 2), which can result in major errors if the input variables are defined inaccurately. Unfortunately, the deviations are expected to be higher at lower I_ta and MATWM (Fig. 2) that are characteristic for permafrost regions. Note that this alternate solution is also implemented in the PERICLIMv1.0 R package.

Appendix C: Validation data

The PERICLIMv1.0 validation was based mostly on previously published data obtained in the period 2011/12 to 2017/18 at four bare permafrost sites located on James Ross Island, northeastern Antarctic Peninsula ($\mathrm{63}{}^{\circ }{\mathrm{49}}^{\prime }$–$\mathrm{63}{}^{\circ }{\mathrm{53}}^{\prime }$ S, $\mathrm{57}{}^{\circ }{\mathrm{50}}^{\prime }$–$\mathrm{57}{}^{\circ }{\mathrm{57}}^{\prime }$ W; 10–340 m a.s.l.; e.g. Hrbáček et al., 2017a, b; Hrbáček and Uxa, 2020), and data collected by the Geophysical Institute Permafrost Laboratory at the University of Alaska Fairbanks in the period 2001/02 to 2016/17 at three vegetated permafrost locations on the coastal plain of the Alaskan Arctic adjacent to the Beaufort Sea ($\mathrm{69}{}^{\circ }{\mathrm{40}}^{\prime }$–$\mathrm{70}{}^{\circ }{\mathrm{22}}^{\prime }$ N, $\mathrm{148}{}^{\circ }{\mathrm{28}}^{\prime }$–$\mathrm{148}{}^{\circ }{\mathrm{43}}^{\prime }$ W; 3–111 m a.s.l.; https://permafrost.gi.alaska.edu/sites_list, last access: 28 June 2019; Romanovsky et al., 2009; Wang et al., 2018) (Table C1). The stations measured air and ground temperatures with thermistor sensors installed in solar radiation shields 1.5 or 2 m above the ground surface and at 6 to 15 depth levels ranging from or near the ground surface to 0.75 m or around 1 m below, and their records were averaged to daily resolution (Romanovsky et al., 2009; Hrbáček et al., 2017a, b; Wang et al., 2018; Hrbáček and Uxa, 2020).

Table C1 Variables used to model air temperature characteristics at the James Ross Island (upper section) and Alaskan Arctic (lower section) validation sites. Note that the superscripts norm and unif indicate normal and uniform distributions describing the respective variables.

See Table 1 for abbreviations.

Active-layer thickness, which corresponds to the maximum annual depth of the 0 ^∘C isotherm (Burn, 1998), was primarily determined by a linear interpolation between the depths of the deepest and the shallowest sensors with maximum annual ground temperatures >0 and ≤0 ^∘C, respectively. Alternatively, it was established by a linear extrapolation of the maximum annual ground temperatures of the two deepest sensors if both were positive. L_t and L_f were defined by a continued prevalence of positive and negative mean daily temperatures, respectively, at the shallowest ground temperature sensor, and for consistency these time windows were also utilized for air temperatures. Since the model assumes that air temperatures are solely positive and negative during the thawing and freezing season, respectively (Fig. 1), positive and negative air temperatures alone were taken to determine MATTS, MATFS, I_ta, and I_fa. Likewise, MAAT was calculated as a length-weighted average of the seasonal air temperature means for a period composed of the thawing season and its preceding freezing season, which is thought to be more representative for active-layer formation than a fixed calendar period (Hrbáček and Uxa, 2020). Annual air temperature range was defined by an annual spread of a 31 d simple central moving average of mean daily air temperatures, with its extremes being considered to substitute MATWM and MATCM. Finally, the ground-surface thawing n factor was derived as a ratio of the thawing index at the shallowest ground temperature sensor and the air thawing index. Consequently, for modelling, the active-layer thickness had to be reduced by the depth of the shallowest ground temperature sensor to ensure consistency of the calculations because the model presumes that the ground-surface thawing n factors transfer between air and ground-surface temperatures (see Hrbáček and Uxa, 2020).

Ground physical properties for the James Ross Island sites were determined in situ or from intact samples collected near the temperature monitoring stations at a depth of 0.1–0.3 m during the thawing seasons 2013/14 to 2018/19 (Hrbáček et al., 2017a; Hrbáček and Uxa, 2020), while those for the Alaskan sites were adapted from Zhang (1993) and Romanovsky and Osterkamp (1997), who took samples from a depth of up to about 0.6 m during the thawing season 1991 and then averaged their characteristics over the full active-layer thickness. Volumetric ground moisture content was established by successive wet and dry weighing or through replicate measurements with time-domain reflectometry probes. Similarly, dry ground bulk density was determined using dry weighing (Zhang, 1993; Hrbáček et al., 2017a; Hrbáček and Uxa, 2020). Ground quartz content was estimated on the basis of the proportions of clay, silt, and sand fractions (see Appendix A), which were ascertained through wet sieving and X-ray diffraction (Hrbáček et al., 2017a; Hrbáček and Uxa, 2020) or assessed visually via soil type (Zhang, 1993). All the substrates were considered fine-grained owing to their relatively high clay–silt contents.

Appendix D: Ground-surface and air thawing index in a two-layer ground

If two distinct ground layers can be distinguished within the palaeo-active layer, the Stefan equation for calculating the active-layer thickness in two-layer ground can be applied. It has been proposed in the following form (Nixon and McRoberts, 1973; Kurylyk, 2015):

$$\begin{array}{}\text{(D1)}& \mathit{\xi }=-Z_{\mathrm{1}}{\displaystyle \frac{k_{{\mathrm{t}}_{\mathrm{2}}}}{k_{{\mathrm{t}}_{\mathrm{1}}}}}+Z_{\mathrm{1}}+\sqrt{{\displaystyle \frac{Z_{\mathrm{1}}^{\mathrm{2}}{k}_{{\mathrm{t}}_{\mathrm{2}}}^{\mathrm{2}}}{k_{{\mathrm{t}}_{\mathrm{1}}}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}{k}_{{\mathrm{t}}_{\mathrm{2}}}{I}_{\text{ts}}}{L{\mathit{\varphi }}_{\mathrm{2}}{\mathit{\rho }}_{\mathrm{w}}}}-{\displaystyle \frac{Z_{\mathrm{1}}^{\mathrm{2}}{k}_{{\mathrm{t}}_{\mathrm{2}}}{\mathit{\varphi }}_{\mathrm{1}}}{k_{{\mathrm{t}}_{\mathrm{1}}}{\mathit{\varphi }}_{\mathrm{2}}}}},\end{array}$$

where Z₁ [m] is the thickness of the top sub-layer, the physical properties of which are indicated by the subscript 1, while the bottom sub-layer is denoted by the subscript 2. The ground-surface index can then be simply expressed from Eq. (D1):

$$\begin{array}{}\text{(D2)}& I_{\text{ts}}={\displaystyle \frac{\left[(\mathit{\xi }+Z_{\mathrm{1}}\frac{k_{{\mathrm{t}}_{\mathrm{2}}}}{k_{{\mathrm{t}}_{\mathrm{1}}}}-Z_{\mathrm{1}}{)}^{\mathrm{2}}-\frac{Z_{\mathrm{1}}^{\mathrm{2}}{k}_{{\mathrm{t}}_{\mathrm{2}}}^{\mathrm{2}}}{k_{{\mathrm{t}}_{\mathrm{1}}}^{\mathrm{2}}}+\frac{Z_{\mathrm{1}}^{\mathrm{2}}{k}_{{\mathrm{t}}_{\mathrm{2}}}{\mathit{\varphi }}_{\mathrm{1}}}{k_{{\mathrm{t}}_{\mathrm{1}}}{\mathit{\varphi }}_{\mathrm{2}}}\right]L{\mathit{\varphi }}_{\mathrm{2}}{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{2}{k}_{{\mathrm{t}}_{\mathrm{2}}}}}.\end{array}$$

As in Eqs. (1) and (2), the product of ϕ and ρ_w can be substituted by that of the gravimetric moisture content and dry bulk density of the ground, but note that the fraction on the far right of Eq. (D1) and at the corresponding place in Eq. (D2) is simplified because the density of water in its numerator and denominator is the same. Subsequent procedures to derive the air temperature characteristics are analogous to those for the one-layer solution (Eqs. 3 to 14).