Articles | Volume 14, issue 4
Geosci. Model Dev., 14, 1865–1884, 2021
https://doi.org/10.5194/gmd-14-1865-2021
Geosci. Model Dev., 14, 1865–1884, 2021
https://doi.org/10.5194/gmd-14-1865-2021

Model description paper 07 Apr 2021

Model description paper | 07 Apr 2021

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

PERICLIMv1.0: a model deriving palaeo-air temperatures from thaw depth in past permafrost regions
Tomáš Uxa1, Marek Křížek2, and Filip Hrbáček3 Tomáš Uxa et al.
  • 1Institute of Geophysics, Czech Academy of Sciences, Prague, Czech Republic
  • 2Department of Physical Geography and Geoecology, Faculty of Science, Charles University, Prague, Czech Republic
  • 3Department of Geography, Faculty of Science, Masaryk University, Brno, Czech Republic

Correspondence: Tomáš Uxa (uxa@ig.cas.cz)

Abstract

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.3C. The past mean annual air temperature modelled experimentally for two sites in the Czech Republic hosting relict cryoturbation structures was between -7.0±1.9 and -3.2±1.5C, 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 (Washburn1980; French2017; Ballantyne2018). 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 Harris1994). Such empirical interpretations are, however, problematic because suitable analogues for past periglacial environments are rare, particularly due to substantial latitudinal contrasts in solar insolation (Williams1975; French2017), 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; Ballantyne2018). Climate controls on the development of most periglacial features are thus poorly known, usually implying broad ranges of climate conditions (Washburn1980; Harris1982, 1994; Karte1983; Wayne1983; Ballantyne and Harris1994; Huijzer and Isarin1997; Ballantyne2018), 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 (Ballantyne2018). 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 Harris1994; Ballantyne2018). 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 Barsch1993; French and Thorn2006).

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 (Williams1961). 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 (French2017). Since the freeze–thaw depth is closely coupled with ground and air temperature conditions (e.g. Frauenfeld et al.2004; Åkerman and Johansson2008; Wu and Zhang2010), it retains a valuable palaeo-climate record that can be approximated based on modern air temperature–freeze–thaw depth relations (Williams1975) or retrieved through an inverse solution of the equations calculating the freeze–thaw depth (Maarleveld1976; French2008). Obviously, this idea is not new, but despite its ingenuity, simplicity, and general acceptance in benchmark periglacial literature (e.g. Washburn1979; Ballantyne and Harris1994; French2017; Ballantyne2018), 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 Williams1975; Maarleveld1976) 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 Nelson2002; Hrbáček and Uxa2020). 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 (Kurylyk2015). 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 Osterkamp1997; Kurylyk and Hayashi2016), but its accuracy is still reasonable (e.g. Klene et al.2001; Shiklomanov and Nelson2002; Hrbáček and Uxa2020). 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 1List of model input and output variables, their symbols, value ranges, and units.

Download Print Version | Download XLSX

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−3], the ground quartz content q [], the ground grain-size class f/c [], the ground-surface thawing n factor nt [], and the annual air temperature range Aa [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 (Lunardini1981):

(1) ξ = 2 k t I ts L ϕ ρ w ,

where kt [Wm-1K-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), Its [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−1] is the specific latent heat of fusion of water, and ρw [1000 kg m−3] is the density of water. Note that Its must be multiplied by the scaling factor of 86 400 s d−1 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.

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f01

Figure 1An idealized course of air temperature during the year described by a sine function with a mean of −4C 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.

Download

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

(2) I ts = ξ 2 L ϕ ρ w 2 k t .

The ground-surface thawing index can then be converted into the air thawing index Ita [C d] through the ground-surface thawing n factor (Lunardini1978), 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):

(3) I ta = I ts n t .

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 Ta(t) [C] can be described by a sine wave (Fig. 1) such as

(4) T a ( t ) = MAAT + A a 2 sin 2 π t P ,

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:

(5) I ta = t 1 t 2 T a ( t ) d t ,

with

(6)t1=arcsin-MAATAa2P2π,(7)t2=π-arcsin-MAATAa2P2π,

where t1 is the time when the air temperature curve crosses the 0 C level from below (∼thawing season begins), while t2 is the time when it crosses this level from above ( thawing season ends) (e.g. Nelson and Outcalt1987). 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 -Aa2<MAAT0. 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 Riseborough2002). 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.

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f02

Figure 2A 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).

Download

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

(8) I fa = MAAT P - I ta .

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

(9)MATWM=MAAT+Aa2,(10)MATCM=MAAT-Aa2.

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

(11)MATTS=ItaLt,(12)MATFS=IfaLf,

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

(13)Lt=t2-t1=π-2arcsin-MAATAa2P2π,(14)Lf=P-Lt.

Unconventionally, the model can also be solved with the range of annual air temperature oscillations defined by MATWM (see Williams1975) if its difference from MAAT (i.e. MATWM−MAAT) is substituted for Aa2 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).

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f03

Figure 3Means 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).

Download

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 ≤2C 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.9C, respectively, and average mean absolute errors achieved 1.1 C for both characteristics (Fig. 3). Their mean absolute deviations were ≤1C at 43 % of the validation sites and ≤2C at 100 % of them. Mean absolute MATWM and MATCM errors reached a maximum of 1.9 and 1.8 C, respectively.

Ita and Ifa were modelled with average mean errors of −60 and −403C 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). Ita was biased by ≤40C d at 29 % of the validation sites and by ≤80C d at 43 % of them. Ifa, which reaches 1 order of magnitude larger values, deviated by ≤400C d at 57 % of the validation sites and by ≤800C d at 100 % of them. Maximum mean absolute departures of Ita and Ifa 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 ≤1C at 100 % of the validation sites, with a maximum mean absolute error of 1.0 C. MATFS deviated by ≤1C at 71 % of the validation sites, and at worst it was 2.5 C.

Since Lt and Lf inherently counteract, their characteristics mirror each other if the values are not rounded. Lt tended to be underestimated by an average of 21 d, while Lf was overestimated by an average of 20 d (Fig. 3), which comprised −20 % and 8 % of the observed average mean Lt and Lf, 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 <-8 and <-4C have been suggested for their formation within coarse- and fine-grained substrates, respectively (Vandenberghe2013; French2017), 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.

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f04

Figure 4Cross 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.

Download

The Brno–Černovice site (491043′′ N, 163856′′ 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 ∼40m 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; Musil1997; 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 (500713′′ N, 122709′′ 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 ∼10m (Š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 (Czudek2005). 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 (Kasse1993; Huijzer and Vandenberghe1998; 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.13m at the Brno–Černovice and Nebanice sites, respectively (Table 2).

Table 2Variables 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.

Download Print Version | Download XLSX

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 cm3 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 (Vandenberghe2013; French2017), we further skewed the moisture content values leftwards using a beta distribution as follows:

(15) ϕ ( x ) = ϕ 0 + f ( x ; α , β ) n - ϕ 0 ,

where ϕ0 [] is the current average volumetric ground moisture content, f(x;α,β) expresses the probability density function of the beta distribution for 0x1 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 Ladanyi2004) 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 (Uxa2017), 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.

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f05

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

Download

Other ground physical properties were supposed to be unchanged since the cryoturbations developed. The dry ground bulk density was 1635±120 and 1645±116kg m−3 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 Johansen1977) 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 Humlum2007) and British Columbia–Yukon and Labrador, Canada (Lewkowicz et al.2012; Way and Lewkowicz2018), 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-Borovskiy1993). 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.6C, 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.6C, respectively, which is within the range of 28–36 C suggested by previous regional estimates for the end of the LGM (Huijzer and Vandenberghe1998) 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 (Carnell2020), which discretizes the probability density functions of the individual model input variables (Table 2) into N non-overlapping intervals of equal probability 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).

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f06

Figure 6Means 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.

Download

MAAT was modelled at -6.6±2.7 to -3.3±2.0C and -7.0±1.9 to -3.2±1.5C at the Brno–Černovice and Nebanice sites, respectively (Fig. 6), which corresponds to its average decline of −16.0 to −12.7C and −15.1 to −11.3C, 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.6C and 4.7±0.8 to 5.5±1.0C at the Brno–Černovice and Nebanice sites, respectively, and culminated with MATWM of 8.4±1.9 to 10.1±2.6C and 7.3±1.3 to 8.5±1.6C (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 %. Ita showed even more consistent outputs as it reached 823±271 to 915±353C d and 704±181 to 721±203C 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, Lt 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 -14.2±1.8 to -9.3±1.5C and -13.7±1.6 to -8.5±1.5C at the Brno–Černovice and Nebanice sites, respectively, and MATCM as low as -23.2±3.2 to -15.0±2.6C and -22.5±2.8 to -13.7±2.5C (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. Ifa spanned -3309±667 to -2041±492C d and -3270±523 to -1873±429C 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, Lf 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 Johansson2008; Wu and Zhang2010). 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 Uxa2020). By contrast, the Alaskan profiles have a two-layer composition with peat over mineral soil (Zhang1993). 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 Osterkamp1997; Kurylyk and Hayashi2016). 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 Hayashi2016). 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 ≤1C, Ita and Ifa deviating on average by −1 % and −8 %, respectively, and Lt and Lf 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 Humlum2007; Lewkowicz et al.2012; Way and Lewkowicz2018) 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 -7.0±1.9 and -3.2±1.5C and its corresponding reduction between −16.0 and −11.3C 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 −12C (Poser1948; Büdel1953; Kaiser1960; Frenzel1967; Goździk1973; Huijzer and Vandenberghe1998; Marks et al.2016). By contrast, it disagrees with slightly milder MAAT reductions of at least −7 and −13 to −6C derived from groundwater data (Corcho Alvarado et al.2011) and borehole temperature logs (Šafanda and Rajver2001), 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 230′′ 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.4C at the study sites, which corresponds to its reduction of between −12.6 and −9C 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.6C 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 Vandenberghe1998; 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 -23.2±3.2 to -13.7±2.5C, which, however, also pertains to other proxy records that show somewhat lower values of −27 to −16.5C (Huijzer and Vandenberghe1998; Marks et al.2016). The GCM-based MATCM exhibits less variability of −21.6 to −17.7C, 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 -7.0±1.9 and -3.2±1.5C is relatively consistent with the MAAT threshold of <-4C commonly suggested for the formation of cryoturbation structures within fine-grained substrates (Vandenberghe2013). At the study sites, however, the cryoturbations consist of coarse-grained materials, for which a MAAT threshold as low as <-8C has been proposed (Vandenberghe2013). This study thus raises questions about the validity of the previously suggested MAAT thresholds for cryoturbation structures (see Vandenberghe2013; French2017) 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).

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f07

Figure 7Sensitivity 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.

Download

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 Harris1994; Matsuoka2011; Ballantyne2018). 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. Karte1983; Matsuoka2001; Ballantyne2018) 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 Harris1994; Ballantyne2018). 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 (Ballantyne1998) 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 Hauck2017).

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.3C. Ita and Ifa deviate on average by −16 % and −12 %, respectively, while Lt and Lf 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 -7.0±1.9 and -3.2±1.5C, and its corresponding reduction was between −16.0 and −11.3C 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 kdry [Wm-1K-1] and saturated ksat [Wm-1K-1] thermal conductivity of a material as follows:

(A1) k t = k dry + ( k sat - k dry ) K e ,

where Ke [] is the dimensionless Kersten number used to normalize the contrast between ksat and kdry based on the degree of saturation.

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

(A2) k dry = 0.135 ρ + 64.7 2700 - 0.947 ρ ,

where the constant 2700 kg m−3 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 (Johansen1977):

(A3) k sat = k s 1 - n k w n ,

where ks [Wm-1K-1] and kw [Wm-1K-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:

(A4) n = 1 - ρ 2700 .

The thermal conductivity of water is set at 0.57 Wm-1K-1, while that of solids is computed as

(A5) k s = k q q k o 1 - q ,

where kq [Wm-1K-1] and ko [Wm-1K-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 Wm-1K-1, while for other minerals it is as follows (Johansen1977):

(A6) k o = 3 , for q < 20 % coarse-grained 2 , otherwise .

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

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 % (Johansen1977). Alternatively, it can also be drawn from a nomogram using the ground texture (Fig. A1).

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f08

Figure A1A 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.

Download

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 (Johansen1977)

(A7) K e = log S + 1 , for S > 10 % fine-grained 0.7 log S + 1 , for S > 5 % coarse-grained ,

where

(A8) S = ϕ n .

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).

https://gmd.copernicus.org/articles/14/1865/2021/gmd-14-1865-2021-f09

Figure A2Thermal 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.

Download

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 Williams1975) if its difference from MAAT (i.e. MATWM−MAAT) is substituted for Aa2 in Eq. (4) and elsewhere. Likewise, MAAT is then calculated numerically using the bisection root-finding algorithm, but it is searched for -<MAAT0 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 Ita 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 Ita 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 (63496353 S, 57505757 W; 10–340 m a.s.l.; e.g. Hrbáček et al.2017a, b; Hrbáček and Uxa2020), 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 (69407022 N, 1482814843 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 Uxa2020).

Table C1Variables 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.

Download Print Version | Download XLSX

Active-layer thickness, which corresponds to the maximum annual depth of the 0 C isotherm (Burn1998), was primarily determined by a linear interpolation between the depths of the deepest and the shallowest sensors with maximum annual ground temperatures >0 and ≤0C, respectively. Alternatively, it was established by a linear extrapolation of the maximum annual ground temperatures of the two deepest sensors if both were positive. Lt and Lf 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, Ita, and Ifa. 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 Uxa2020). 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 Uxa2020).

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 Uxa2020), 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 (Zhang1993; Hrbáček et al.2017a; Hrbáček and Uxa2020). 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 Uxa2020) or assessed visually via soil type (Zhang1993). 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 McRoberts1973; Kurylyk2015):

(D1) ξ = - Z 1 k t 2 k t 1 + Z 1 + Z 1 2 k t 2 2 k t 1 2 + 2 k t 2 I ts L ϕ 2 ρ w - Z 1 2 k t 2 ϕ 1 k t 1 ϕ 2 ,

where Z1 [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):

(D2) I ts = ( ξ + Z 1 k t 2 k t 1 - Z 1 ) 2 - Z 1 2 k t 2 2 k t 1 2 + Z 1 2 k t 2 ϕ 1 k t 1 ϕ 2 L ϕ 2 ρ w 2 k t 2 .

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).

Code and data availability

The latest version of PERICLIMv1.0 is available as an R package from https://github.com/tomasuxa/PERICLIMv1.0 (last access: 15 January 2021) under the GPLv3 license. The exact version of the model used to produce this paper is archived at https://doi.org/10.5281/zenodo.4562435 (Uxa2021). The validation datasets from James Ross Island are available upon request from Filip Hrbáček (hrbacekfilip@gmail.com), whereas those from the Alaskan Arctic can be retrieved from https://permafrost.gi.alaska.edu/sites_list (last access: 28 June 2019).

Author contributions

TU came up with an initial idea with feedback from MK, developed the model, and evaluated it against data from James Ross Island and the Alaskan Arctic, which were processed by FH and TU, respectively. TU then tested it for the derivation of palaeo-air temperature characteristics using data sampled collectively with MK in the Czech Republic. TU drew figures and wrote the paper with inputs from MK and FH. All authors reviewed and approved the final version of the paper.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We thank Tereza Dlabáčková for her assistance at the Nebanice site and with sample analysis. The Geophysical Institute Permafrost Laboratory at the University of Alaska Fairbanks is acknowledged for its continuous effort in collecting temperature data across Alaska and their online dissemination. Last but not least, we thank the reviewers and the editor, who significantly contributed to the improvement of the paper.

Financial support

The PERICLIMv1.0 development and evaluation were supported by the Czech Science Foundation, project number 17-21612S. The validation datasets from James Ross Island were collected thanks to the Ministry of Education, Youth and Sports of the Czech Republic, project number VAN2020/01.

Review statement

This paper was edited by Andrew Wickert and reviewed by four anonymous referees.

References

Åkerman, H. J. and Johansson, M.: Thawing permafrost and thicker active layers in sub-arctic Sweden, Permafrost Periglac., 19, 279–292, https://doi.org/10.1002/ppp.626, 2008. a, b

Andersland, O. B. and Ladanyi, B.: Frozen Ground Engineering, 2nd Edition, John Wiley & Sons, Hoboken, USA, 2004. a

Andrieux, E., Bateman, M. D., and Bertran, P.: The chronology of Late Pleistocene thermal contraction cracking derived from sand wedge OSL dating in central and southern France, Global Planet. Change, 162, 84–100, https://doi.org/10.1016/j.gloplacha.2018.01.012, 2018. a

Balatka, B., Kalvoda, J., Steklá, T., and Štěpančíková, P.: Morphostratigraphy of river terraces in the Eger valley (Czechia) focused on the Smrčiny Mountains, the Chebská pánev Basin and the Sokolovská pánev Basin, AUC Geogr., 54, 240–259. https://doi.org/10.14712/23361980.2019.21, 2019. a, b

Ballantyne, C. K.: Age and Significance of Mountain-Top Detritus, Permafrost Periglac., 9, 327–345, https://doi.org/10.1002/(SICI)1099-1530(199810/12)9:4<327::AID-PPP298>3.0.CO;2-9, 1998. a

Ballantyne, C. K.: Periglacial Geomorphology, John Wiley & Sons, Hoboken, USA, 2018. a, b, c, d, e, f, g, h, i

Ballantyne, C. K. and Harris, C.: The Periglaciation of Great Britain, Cambridge University Press, Cambridge, UK, 1994. a, b, c, d, e, f

Barsch, D.: Periglacial geomorphology in the 21st century, Geomorphology, 7, 141–163, https://doi.org/10.1016/B978-0-444-89971-2.50011-0, 1993. a

Bertran, P., Andrieux, E., Antoine, P., Coutard, S., Deschodt, L., Gardère, P., Hernandez, M., Legentil, C., Lenoble, A., Liard, M., Mercier, N., Moine, O., Sitzia, L., and Van Vliet-Lanoë, B.: Distribution and chronology of Pleistocene permafrost features in France: database and first results, Boreas, 43, 699–711, https://doi.org/10.1111/bor.12025, 2014. a

Bubík, M., Hanžl, P., Havlíček, P., Novák, Z., Otava, J., Petrová, J., Valoch, K., and Vít, J.: Výběr některých zajímavých lokalit – 69 Černovice (B4) [Selection of some interesting locations – 69 Černovice (B4)], in: Geologie Brna a okolí [Geology of Brno and surroundings], edited by: Müller, P. and Novák, Z., Český geologický ústav, Praha, Czech Republic, 83, 2000. a

Büdel, J.: Die “periglazial”-morphologisehen Wirkungen des Eiszeitklimas auf der ganzen Erde, Erdkunde, 7, 249–256, 1953. a

Burn, C. R.: The Active Layer: Two Contrasting Definitions, Permafrost Periglac., 9, 411–416, https://doi.org/10.1002/(SICI)1099-1530(199810/12)9:4<411::AID-PPP292>3.0.CO;2-6, 1998. a

Carnell, R.: lhs: Latin Hypercube Samples, R package version 1.0.2, https://cran.r-project.org/web/packages/lhs, last access: 8 August 2020. a

Corcho Alvarado, J. A., Leuenberger, M., Kipfer, R., Paces, T., and Purtschert, R.: Reconstruction of past climate conditions over central Europe from groundwater data, Quaternary Sci. Rev., 30, 3423–3429, https://doi.org/10.1016/j.quascirev.2011.09.003, 2011. a

Czudek, T.: Vývoj reliéfu krajiny České republiky v kvartéru [Quaternary Development of Landscape Relief of the Czech Republic], Moravské zemské muzeum, Brno, Czech Republic, 2005. a

Dong, Y., McCartney, J. S., and Lu, N.: Critical Review of Thermal Conductivity Models for Unsaturated Soils, Geotech. Geol. Eng., 33, 207–221. https://doi.org/10.1007/s10706-015-9843-2, 2015. a, b

Engel, Z., Křížek, M., Braucher, R., Uxa, T., Krause, D., and AsterTeam: 10Be exposure age for sorted polygons in the Sudetes Mountains, Permafrost Periglac., 32, 154–168, https://doi.org/10.1002/ppp.2091, 2021. a

Frauenfeld, O. W., Zhang, T., Barry, R. G., and Gilichinsky, D.: Interdecadal changes in seasonal freeze and thaw depths in Russia, J. Geophys. Res.-Atmos., 109, D05101, https://doi.org/10.1029/2003JD004245, 2004. a, b

French, H.: Recent Contributions to the Study of Past Permafrost, Permafrost Periglac., 19, 179–194, https://doi.org/10.1002/ppp.614, 2008. a

French, H. and Thorn, C. E.: The changing nature of periglacial geomorphology, Geomorphologie, 12, 165–174, https://doi.org/10.4000/geomorphologie.119, 2006. a

French, H. M.: The Periglacial Environment, 4th Edition, John Wiley & Sons, Hoboken, USA, 2017. a, b, c, d, e, f, g

Frenzel, B.: Die Klimaschwankungen des Eiszeitalters. Friedr. Vieweg & Sohn, Braunschweig, Germany, 1967. a

Gisnås, K., Westermann, S., Schuler, T. V., Melvold, K., and Etzelmüller, B.: Small-scale variation of snow in a regional permafrost model, The Cryosphere, 10, 1201–1215, https://doi.org/10.5194/tc-10-1201-2016, 2016. a

Gisnås, K., Etzelmüller, B., Lussana, C., Hjort, J., Sannel A. B. K., Isaksen, K., Westermann, S., Kuhry, P., Christiansen, H. H., Frampton, A., and Åkerman, J.: Permafrost Map for Norway, Sweden and Finland, Permafrost Periglac., 28, 359–378, https://doi.org/10.1002/ppp.1922, 2017. a

Goździk, J.: Geneza i pozycja stratygraficzna struktur peryglacjalnych w środkowej Polsce [The genesis and stratigraphical position of periglacial structures in central Poland], Acta Geogr. Lodz., 31, 5–117, 1973. a

Harris, S. A.: Identification of permafrost zones using selected permafrost landforms, in: Proceedings of the 4th Canadian Permafrost Conference, Calgary, Canada, 2–6 March 1981, 49–58, 1982. a

Harris, S. A.: Climatic Zonality of Periglacial Landforms in Mountain Areas, Arctic, 47, 184–192, https://doi.org/10.14430/arctic1288, 1994. a

He, H., Zhao, Y., Dyck, M. F., Si, B., Jin, H., Lv, J., and Wang, J.: A modified normalized model for predicting effective soil thermal conductivity, Acta Geotech., 12, 1281–1300. https://doi.org/10.1007/s11440-017-0563-z, 2017. a, b

Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., and Jarvis, A.: Very high resolution interpolated climate surfaces for global land areas, Int. J. Climatol., 25, 1965–1978, https://doi.org/10.1002/joc.1276, 2005. a

Hrbáček, F. and Uxa, T.: The evolution of a near-surface ground thermal regime and modeled active-layer thickness on James Ross Island, Eastern Antarctic Peninsula, in 2006–2016, Permafrost Periglac., 31, 141–155, https://doi.org/10.1002/ppp.2018, 2020. a, b, c, d, e, f, g, h, i, j

Hrbáček, F., Kňažková, M., Nývlt, D., Láska, K., Mueller, C. W., and Ondruch, J.: Active layer monitoring at CALM-S site near J. G. Mendel Station, James Ross Island, eastern Antarctic Peninsula, Sci. Total Environ., 601–602, 987–997, https://doi.org/10.1016/j.scitotenv.2017.05.266, 2017a. a, b, c, d, e, f

Hrbáček, F., Nývlt, D., and Láska, K.: Active layer thermal dynamics at two lithologically different sites on James Ross Island, Eastern Antarctic Peninsula, Catena, 149, 592–602, https://doi.org/10.1016/j.catena.2016.06.020, 2017b. a, b

Huijzer, A. S. and Isarin, R. F. B.: The reconstruction of past climates using multi-proxy evidence: An example of the weichselian pleniglacial in northwest and central Europe, Quaternary Sci. Rev., 16, 513–533, https://doi.org/10.1016/S0277-3791(96)00080-7, 1997. a

Huijzer, B. and Vandenberghe, J.: Climatic reconstruction of the Weichselian Pleniglacial in northwestern and central Europe, J. Quaternary Sci., 13, 391–417, https://doi.org/10.1002/(SICI)1099-1417(1998090)13:5<391::AID-JQS397>3.0.CO;2-6, 1998. a, b, c, d, e

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, https://doi.org/10.5194/tc-6-613-2012, 2012. a

Johansen, Ø.: Thermal conductivity of soils, Draft Translation 637, US Army Cold Regions Research and Engineering Laboratory, Hanover, USA, 1977. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q

Juliussen, H. and Humlum, O.: Towards a TTOP Ground Temperature Model for Mountainous Terrain in Central-Eastern Norway, Permafrost Periglac., 18, 161–184, https://doi.org/10.1002/ppp.586, 2007. a, b

Kaiser, K.: Klimazeugen des periglazialen Dauerfrostbodens in Mittel- und Westeuropa. Ein Beitrag zur Rekonstruktion des Klimas der Glaziale des quartären Eiszeitalters, E & G Quaternary Sci. J., 11, 121–141, https://doi.org/10.23689/fidgeo-1214, 1960. a

Karte, J.: Periglacial Phenomena and their Significance as Climatic and Edaphic Indicators, Geoj., 7, 329–340, https://doi.org/10.1007/BF00241455, 1983. a, b

Kasse, C.: Periglacial environments and climatic development during the Early Pleistocene Tiglian stage (Beerse Glacial) in northern Belgium, Geologie en Mijnbouw, 72, 107–123, 1993. a

Kasse, C., Vandenberghe, J., Van Huissteden, J., Bohncke, S. J. P., and Bos, J. A. A.: Sensitivity of Weichselian fluvial systems to climate change (Nochten mine, eastern Germany), Quaternary Sci. Rev., 22, 2141–2156, https://doi.org/10.1016/S0277-3791(03)00146-X, 2003. a

Klene, A. E., Nelson, F. E., Shiklomanov, N. I., and Hinkel, K. M.: The N-Factor in Natural Landscapes: Variability of Air and Soil-Surface Temperatures, Kuparuk River Basin, Alaska, USA, Arct. Antarct. Alp. Res., 33, 140–148, https://doi.org/10.2307/1552214, 2001. a, b, c

Kudryavtsev, V. A., Garagulya, L. S., Kondratyeva, K. A., and Melamed, V. G.: Fundamentals of Frost Forecasting in Geological Engineering Investigations, Draft Translation 606, US Army Cold Regions Research and Engineering Laboratory, Hanover, USA, 1977. a

Kuneš, P., Pelánková, B., Chytrý, M., Jankovská, V., Pokorný, P., and Petr, L.: Interpretation of the last-glacial vegetation of eastern-central Europe using modern analogues from southern Siberia, J. Biogeogr., 35, 2223–2236, https://doi.org/10.1111/j.1365-2699.2008.01974.x, 2008. a

Kurylyk, B. L.: Discussion of “A Simple Thaw-Freeze Algorithm for a Multi-Layered Soil using the Stefan Equation” by Xie and Gough (2013), Permafrost Periglac., 26, 200–206, https://doi.org/10.1002/ppp.1834, 2015. a, b

Kurylyk, B. L. and Hayashi, M.: Improved Stefan Equation Correction Factors to Accommodate Sensible Heat Storage during Soil Freezing or Thawing, Permafrost Periglac., 27, 189–203, https://doi.org/10.1002/ppp.1865, 2016. a, b, c

Lewkowicz, A. G., Bonnaventure, P. P., Smith, S. L., and Kuntz, Z.: Spatial and thermal characteristics of mountain permafrost, Northwest Canada, Geogr. Ann. A, 94, 195–213, https://doi.org/10.1111/j.1468-0459.2012.00462.x, 2012. a, b

Lindgren, A., Hugelius, G., Kuhry, P., Christensen, T. R., and Vandenberghe, J.: GIS-based Maps and Area Estimates of Northern Hemisphere Permafrost Extent during the Last Glacial Maximum, Permafrost Periglac., 27, 6–16, https://doi.org/10.1002/ppp.1851, 2016. a

Lunardini, V. J.: Theory of N-factors and correlation of data, in: Proceedings of the 3rd International Conference on Permafrost, Vol. 1, Edmonton, Canada, 10–13 July 1978, 40–46, 1978. a

Lunardini, V. J.: Heat Transfer in Cold Climates, Van Nostrand Reinhold Co., New York, USA, 1981. a

Maarleveld, G. C.: Periglacial phenomena and the mean annual temperature during the last glacial time in the Netherlands, Biul. Peryglac., 26, 57–78, 1976. a, b

Marks, L., Gałązka, D., and Woronko, B.: Climate, environment and stratigraphy of the last Pleistocene glacial stage in Poland, Quatern. Int., 420, 259–271, https://doi.org/10.1016/j.quaint.2015.07.047, 2016. a, b, c

Matsuoka, N.: Solifluction rates, processes and landforms: a global review, Earth-Sci. Rev., 55, 107–134, https://doi.org/10.1016/S0012-8252(01)00057-5, 2001. a

Matsuoka, N.: Climate and material controls on periglacial soil processes: Toward improving periglacial climate indicators, Quaternary Res., 75, 356–365, https://doi.org/10.1016/j.yqres.2010.12.014, 2011. a

McKay, M. D., Beckman, R. J., and Conover, W. J.: Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code, Technometrics, 21, 239–245, https://doi.org/10.1080/00401706.1979.10489755, 1979. a, b

Musil, R.: Tuřanská terasa Svitavy v Brně [Fluvial Tuřany terrace of the Svitava River in Brno], Geologické výzkumy na Moravě a ve Slezsku, 4, 14–17, 1997. a

Musil, R., Karásek, J., Seitl, L., and Valoch, K.: Fluviální akumulace v Brně–Černovicích [Fluvial aggradational terraces at Brno–Černovice], Geologické výzkumy na Moravě a ve Slezsku, 3, 28–31, 1996. a

Nelson, F. E. and Outcalt, S. I.: A Computational Method for Prediction and Regionalization of Permafrost, Arctic Alpine Res., 19, 279–288, https://doi.org/10.2307/1551363, 1987. a

Nixon, J. F. and McRoberts, E. C.: A Study of Some Factors Affecting the Thawing of Frozen Soils, Can. Geotech. J., 10, 439–452, https://doi.org/10.1139/t73-037, 1973. a

Nyland, K. E., Nelson, F. E., Figueiredo, P. M.: Cosmogenic 10Be and 36Cl geochronology of cryoplanation terraces in the Alaskan Yukon-Tanana Upland, Quaternary Res., 97, 157–166, https://doi.org/10.1017/qua.2020.25, 2020. a

Poser, H.: Boden- und Klimaverhältnisse in Mittel- und Westeuropa während der Würmeiszeit, Erdkunde, 2, 53–68, 1948. 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, https://doi.org/10.1002/(SICI)1099-1530(199701)8:1%3C1::AID-PPP243%3E3.0.CO;2-U, 1997. a, b, c

Romanovsky, V. E., Kholodov, A. L., Cable, W. L., Cohen, L., Panda, S., Marchenko, S., Muskett, R. R., and Nicolsky, D.: Network of Permafrost Observatories in North America and Russia, NSF Arctic Data Center, Santa Barbara, CA, USA, https://doi.org/10.18739/A2SH27, 2009. a, b

Šafanda, J. and Rajver, D.: Signature of the last ice age in the present subsurface temperatures in the Czech Republic and Slovenia, Glob. Planet. Change, 29, 241–257, https://doi.org/10.1016/S0921-8181(01)00093-5, 2001. a

Šantrůček, P., Králík, F., and Kvičinský, Z.: Geologická mapa ČR 1 : 50 000, List 11–14 Cheb [Geological map of the Czech Republic 1:50 000, Sheet 11–14 Cheb], Český geologický ústav, Praha, Czech Republic, 1994. a

Shiklomanov, N. I. and Nelson, F. E.: Active-Layer Mapping at Regional Scales: A 13-Year Spatial Time Series for the Kuparuk Region, North-Central Alaska, Permafrost Periglac., 13, 219–230, https://doi.org/10.1002/ppp.425, 2002. a, b

Shur, Y., Hinkel, K. M., and Nelson, F. E.: The Transient Layer: Implications for Geocryology and Climate-Change Science. Permafrost Periglac., 16, 5–17, https://doi.org/10.1002/ppp.518, 2005. a

Shur, Y. L. and Slavin-Borovskiy, V. B.: N-factor maps of Russian permafrost region, in: Proceedings of the 6th International Conference on Permafrost, Vol. 1, Beijing, China, 5–9 July 1993, 564–568, 1993. a

Smith, M. W. and Riseborough, D. W.: Climate and the Limits of Permafrost: A Zonal Analysis, Permafrost Periglac., 13, 1–15, https://doi.org/10.1002/ppp.410, 2002. a

Špičáková, L., Uličný, D., and Koudelková, G.: Tectonosedimentary Evolution of the Cheb basin (NW Bohemia, Czech Republic) between Late Oligocene and Pliocene: a preliminary note, Stud. Geophys. Geod., 44, 556–580, https://doi.org/10.1023/A:1021819802569, 2000. a

Stefan J.: Über die Theorie der Eisbildung, insbesondere über die Eisbildung im Polarmeere, Ann. Phys., 278, 269–286, https://doi.org/10.1002/andp.18912780206, 1891. a

Uxa, T.: Discussion on “Active Layer Thickness Prediction on the Western Antarctic Peninsula” by Wilhelm et al. (2015), Permafrost Periglac., 28, 493–498, https://doi.org/10.1002/ppp.1888, 2017. a

Uxa, T.: PERICLIMv1.0: A model deriving palaeo-air temperatures from thaw depth in past permafrost regions (Version 1.0), Zenodo, https://doi.org/10.5281/zenodo.4562435, 2021. a

Uxa, T., Mida, P., and Křížek, M.: Effect of Climate on Morphology and Development of Sorted Circles and Polygons, Permafrost Periglac., 28, 663–674, https://doi.org/10.1002/ppp.1949, 2017. a

Vandenberghe, J.: Cryoturbation structures, in: Encyclopedia of Quaternary Science, 2nd Edition, edited by: Elias, S. A. and Mock, C. J., Elsevier, Amsterdam, the Netherlands, 430–435, https://doi.org/10.1016/B978-0-444-53643-3.00096-0, 2013. a, b, c, d, e

Vandenberghe, J., French, H. M., Gorbunov, A., Marchenko, S., Velichko, A. A., Jin, H., Cui, Z., Zhang, T., and Wan, X.: The Last Permafrost Maximum (LPM) map of the Northern Hemisphere: permafrost extent and mean annual air temperatures, 25–17 ka BP, Boreas, 43, 652–666, https://doi.org/10.1111/bor.12070, 2014. a

Wang, K., Jafarov, E., Overeem, I., Romanovsky, V., Schaefer, K., Clow, G., Urban, F., Cable, W., Piper, M., Schwalm, C., Zhang, T., Kholodov, A., Sousanes, P., Loso, M., and Hill, K.: A synthesis dataset of permafrost-affected soil thermal conditions for Alaska, USA, Earth Syst. Sci. Data, 10, 2311–2328, https://doi.org/10.5194/essd-10-2311-2018, 2018. a, b

Washburn, A. L.: Geocryology: A Survey of Periglacial Environments, Edward Arnold, London, UK, 1979. a

Washburn, A. L.: Permafrost features as evidence of climatic change, Earth-Sci. Rev., 15, 327–402, https://doi.org/10.1016/0012-8252(80)90114-2, 1980. a, b

Way, R. G. and Lewkowicz, A. G.: Environmental controls on ground temperature and permafrost in Labrador, northeast Canada, Permafrost Periglac., 29, 73–85, https://doi.org/10.1002/ppp.1972, 2018. a, b

Wayne, W. J.: Paleoclimatic Inferences from Relict Cryogenic Features in Alpine Regions, in: Proceedings of the 4th International Conference on Permafrost, Vol. 1, Fairbanks, USA, 17–22 July 1983, 1378–1383, 1983. a

Westermann, S., Elberling, B., Højlund Pedersen, S., Stendel, M., Hansen, B. U., and Liston, G. E.: Future permafrost conditions along environmental gradients in Zackenberg, Greenland, The Cryosphere, 9, 719–735, https://doi.org/10.5194/tc-9-719-2015, 2015. 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, https://doi.org/10.5194/gmd-9-523-2016, 2016. a

Wicky, J. and Hauck, C.: Numerical modelling of convective heat transport by air flow in permafrost talus slopes, The Cryosphere, 11, 1311–1325, https://doi.org/10.5194/tc-11-1311-2017, 2017.  a

Williams, P. J.: Climatic Factors Controlling the Distribution of Certain Frozen Ground Phenomena, Geogr. Ann., 43, 339–347, https://doi.org/10.1080/20014422.1961.11880994, 1961. a

Williams, R. B. G.: The British climate during the last glaciation: an interpretation based on periglacial phenomena, in: Ice Ages Ancient and Modern, Wright, A. E. and Moseley, F. (Eds.), Seel House, Liverpool, UK, 1975. a, b, c, d, e

Wu, Q. and Zhang, T.: Changes in active layer thickness over the Qinghai-Tibetan Plateau from 1995 to 2007. J. Geophys. Res.-Atmos., 115, D09107, https://doi.org/10.1029/2009JD012974, 2010. a, b

Zhang, M., Bi, J., Chen, W., Zhang, X., and Lu, J.: Evaluation of calculation models for the thermal conductivity of soils, Int. Commun. Heat Mass, 94, 14–23, https://doi.org/10.1016/j.icheatmasstransfer.2018.02.005, 2018. a, b

Zhang, T.: Climate, Seasonal Snow Cover and Permafrost Temperatures in Alaska North of the Brooks Range, Ph.D. thesis, University of Alaska, Fairbanks, USA, 232 pp., 1993. a, b, c, d

Download
Short summary
We present a simple 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. Its evaluation against modern temperature records and an experimental palaeo-air temperature reconstruction showed relatively high model accuracy, which suggests that it could become a useful tool for reconstructing Quaternary palaeo-environments.