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, or blockfields, are among the most common relics of colder climates, which repetitively occurred throughout the Quaternary, and, 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 been far from reliable. This study introduces and evaluates a new simple inverse modelling scheme PERICLIMv1.0 5 (PERIglacial CLIMate) that aims to overcome these flaws through deriving 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. Besides, 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 ◦C and −3.2±1.5 ◦C and its corresponding reduction was 10 between −16.0 ◦C and −11.3 ◦C in comparison with the 1981–2010 period, which is well in line with earlier reconstructions utilizing various palaeo-archives. The relatively high model success rate is promising and suggests that it could become a useful tool for reconstructing Quaternary palaeo-environments across vast areas of midand low-latitudes where relict periglacial assemblages frequently occur, but their full potential remains to be exploited.


15
Many mid-and low-latitude regions of the world host a number of distinctive landforms and subsurface structures collectively termed as relict periglacial features, such as various kinds of patterned ground, cryoturbations, frost wedges, solifluction structures, or 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 rest in places where other palaeo-indicators are rare or absent, or emerged at different times, which enhances their relevance as archives of past environmental conditions. So far, 20 relict periglacial features have been used to reconstruct past climates almost exclusively on the basis of climate thresholds of their active counterparts that are mostly situated in present-day high-latitude periglacial environments (Ballantyne and Harris, of modelling tools and foster the development of new quantitative methods in palaeo-environmental reconstructions utilizing relict periglacial features in order to raise their reputation as palaeo-proxy indicators.

Model description
The PERICLIMv1.0 builds on an inverse solution of the Stefan (1891) equation, which has originally been developed to 60 determine the thickness of sea ice, but it also well describes the thaw propagation in ice-bearing grounds, 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 latent heat is the only energy sink associated with its thawing, that is, 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 as well as 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). Besides, its simplicity and low requirements for input data as compared to more complex analytical or numerical models (e.g., Kudryavtsev et al., 1977;Jafarov et al., 2012;Westermann et al., 2016) is highly advantageous for 70 palaeo-applications as fewer assumptions have to be made. Here, it is solved for a 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 as well as freezing-season air temperature attributes as detailed below. Please note that its code is implemented and disseminated as R package.

Driving parameters
The PERICLIMv1.0 is driven by 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 ('fine' or 'coarse'), 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 MATWM and coldest MATCM month ( • C). The ground physical 80 parameters are assumed to characterize the entire modelling domain (∼active layer) and are treated as constant over time.

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): 85 where k t is the thermal conductivity of the thawed ground (W m −1 K −1 ) calculated here as a function of its dry bulk density and volumetric moisture content using the Johansen's (1977) thermal-conductivity model (Appendix A), I ts is the groundsurface thawing index defined as a sum of positive daily ground-surface temperatures in the thawing season ( • C d) (Fig. 1), L is the specific latent heat of fusion of water (334 000 J kg −1 ), and ρ w is the density of water (1 000 kg m −3 ). Please note that I ts must be multiplied by the scaling factor of 86 400 s d −1 in Eq.
(1) to obtain the active-layer thickness in meters. Besides, 90 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.
The ground-surface thawing index required to reach the specific active-layer thickness can be obtained if Eq. (1) is rearranged such as:  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): Note that the effect of snow cover on ground-surface temperatures does not have to be accounted for because the air thawing 100 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 low. However, it should be borne in mind that it is suitable especially for locations without long-lasting snow cover, which might disrupt the coupling between the air and ground temperatures during the thawing season (e.g., Gisnås et al., 2016).

105
The temporal evolution of air temperature over the year T a (t) ( • C) can be well described by a sine wave ( Fig. 1) such as: where t is the time (d) and P is the period of air temperature oscillations (365 d). 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 (Appendix B).
The air thawing index represents positive area under the annual air temperature curve ( Fig. 1) and can be calculated by integrating Eq. (4) over the thawing season as follows: where t 1 is the time when the air temperature curve crosses the zero-degree Celsius level from below (∼thawing season begins), while t 2 is the time when it crosses this level from above (∼thawing season ends) (e.g., Nelson and Outcalt, 1987).

120
Unfortunately, Eq. (5) has no analytical solution for MAAT. The latter can be derived from a nomogram (Fig. 2), but here it is calculated numerically using the bisection root-finding method searching for MAAT such that −A a < MAAT ≤ 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 permafrost-seasonal frost boundary usually occurs at slightly negative MAAT (Smith and Riseborough,125 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.
Once MAAT is known, the air freezing index I fa ( • C d) can be simply computed as: Furthermore, MATWM and MATCM is calculated as: Mean air temperature of the thawing MATTS and freezing MATFS season ( • C) is defined as: where L t and L f is the duration of the thawing and freezing season, respectively (d), which is expressed from Eq. (6) and (7) as: Unsurprisingly, but importantly, solutions based on alternate driving parameters, as suggested above, produce identical out-140 comes, allowing the model adaptations to specific situations and available data.

Model validation
Performance of the model was tested using the same simulation schemes as for palaeo-applications detailed in the next section (see Sect. 4.2.5), 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 agreements and clear trends 145 along the identity lines 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 outcomes tended to be slightly underestimated, with those from James Ross Island being somewhat more accurate and less scattered than those from Alaska (Fig. 3).
MAAT exhibited an average mean error and an average mean absolute error of −1.3 • C and 1.3 • C, respectively (Fig. 3).
Mean absolute error was ≤ 1 • C and ≤ 2 • C at 57 % and 86 % of the validation sites, respectively, and maximally was 2.4 • C.
MATWM and MATCM showed slightly lower bias as their average mean errors attained −1.0 • C 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 were at worst 1.9 • C and 1.8 • C, respectively.
I ta and I fa were modelled with average mean errors of −60 • C d and −403 • C d, respectively, which corresponds to −16 % 155 and −12 % of the observed average mean values, and their average mean absolute errors reached 71 • C d and 403 • C d, respectively (Fig. 3). I ta biased by ≤ 40 • C d at 29 % of the validation sites and by ≤ 80 • C d at 43 % of them. I fa , which achieves one order of magnitude larger values, deviated by ≤ 400 • C d at 57 % of the validations and by ≤ 800 • C d at 100 % of them.
Maximum mean absolute departures of I ta and I fa achieved 106 • C d and 789 • C d, respectively.
MATTS exhibited an average mean error and an average mean absolute error of −0.5 • C and 0.5 • C, respectively, while for 160 MATFS the mean errors averaged −0.6 • C and 1.1 • C (Fig. 3). The agreement between the modelled and observed MATTS was ≤ 1 • C at 100 % of the validation sites, with 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

Model application
As a feasibility study, the model was utilized experimentally for derivation of palaeo-air temperature conditions on the basis of 170 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 as well as the presence of permafrost at the time of their development. Also, MAAT thresholds of <−8 • C and <−4 • C have been suggested for their 175 formation within coarse-and fine-grained substrates, respectively (Vandenberghe, 2013;French, 2017), the validity of which can be assessed with the model as well.

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 those attributes that allowed to define the most plausible ranges of the model driving parameters.  , which is located about 1.5 km east of the present channelized river bed and ∼40 m higher. The material tends to coarsen 185 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 (50 • 07 13 N, 12 • 27 09 E, 430 m asl) is a currently inactive sand-gravel pit about 1.5 km west of the 190 centre of Nebanice, West 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 (∼MIS 10), and now is situated about 200 m north of the river bed 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, 195 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 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) and arose at the ends of the major cold events when permafrost started to degrade (sensu Vandenberghe, 2013). We believe that the studied relict cryoturbations indicate past environmental conditions similar to those 200 at the end of the Last Glacial Maximum (LGM) when analogous features formed in the nearby regions (e.g., Kasse, 1993;Kasse et al., 2003;Bertran et al., 2014). Table 2. Driving parameters used to model palaeo-air temperature characteristics related to the cryoturbations at the Brno-Černovice and Nebanice site. Note that the parameters described by means and standard deviations (ξ, ρ, nt, Aa) had normal distributions, while those defined by simple ranges were assumed to be beta-(φ) or uniformly distributed (q). The palaeo-active-layer thickness was considered to be represented by the vertical extent of the cryoturbated horizons ( Fig. 4), 205 which was determined in horizontal steps of 0.2 m along the entire length of each cross-section (4.0-8.2 m), and attained 1.58±0.28 and 1.40±0.13 m at the Brno-Černovice and Nebanice site, respectively (Table 2).

Ground physical properties
Ground physical properties were chosen using intact samples collected at four representative positions along four vertical profiles within each cross-section (16 samples per cross-section) into 100 cm 3 stainless steel cylinders, which were then weighted 210 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 site, respectively, but because cryoturbations require well saturated conditions for their development (Vandenberghe, 2013;French, 2017), these values were assumed to be the lower moisture thresholds, while 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 achieved 39.4 % and 39.1 %, respectively.

215
Moisture, however, usually tends to be close to saturation for at least part of the thawing season (Woo, 2012), and thus we further skewed its values leftwards using a beta distribution as follows: where φ 0 is the current average volumetric ground moisture content (−), f (x; α, β) expresses the probability density function of the beta distribution for 0 ≤ x ≤ 1 with shape parameters tentatively assumed to be α = 5 and β = 2 that yield a mode 220 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 site, respectively, which corresponds to the 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 an unfrozen moisture. Slightly lower moisture content of the eastern site in Brno-Černovice as 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.
Contrarily, the other ground physical properties were supposed to be representative for former conditions as such. The dry 230 ground bulk density was 1635±120 kg m −3 and 1645±116 kg m −3 at the Brno-Černovice and Nebanice site, 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).

Ground-surface thawing n-factor 235
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 little vegetated (sparse lichens, mosses, or grasses) surfaces, which are assumed to have existed at the study sites when the cryoturbations originated because treeless landscapes dominated in South Moravia and West Bohemia at that time (e.g., Kuneš et al., 2008). We took over a total of forty-one ground-surface 240 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 as 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 245 1.03±0.12, which was assigned to both study sites (Table 2).

Annual air temperature range
As a starting point, we used the annual air temperature ranges based on monthly air temperatures measured  (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).

255
The model was solved by the Monte Carlo simulation with Latin hypercube sampling (LHS) (McKay et al., 1979) using lhs R package (Cannel, 2020), which discretizes the probability density functions of the individual driving parameters (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 260 driving parameters, and thus it requires fewer simulations to achieve stable outputs if N is large enough.
We realized 1 000 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 driving parameters and provided the most plausible palaeo-air temperature characteristics that account for the uncertainty as well as natural temporal variability of the inputs.

Modelled palaeo-air temperature characteristics 265
Most modelled palaeo-air temperature characteristics changed in 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 scatters of 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).
MAAT was modelled at −6.6±2.7 to −3.3±2.0 • C and −7.0±1.9 to −3.2±1.5 • C at the Brno-Černovice and Nebanice site, 270 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, respectively, were 3.3 • C and 3.8 • C at the Brno-Černovice and Nebanice site, 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 site, 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 attained 823±271 to 915±353 • C d and 704±181 to 721±203 • C d at the Brno-Černovice and Nebanice site, respectively (Fig. 6), and, as such,
Obviously, the freezing seasons exhibited the largest scatters of the model outputs as well as the highest variability among the scenarios with MATFS of −14.2±1.8 to −9.3±1.5 • C and −13.7±1.6 to −8.5±1.5 • C at the Brno-Černovice and Nebanice site, respectively, and MATCM as low as −23.2±3.2 to −15.0±2.6 • C and −22.5±2.8 to −13.7±2.5 • C (Fig. 6). It thus follows 285 that their variations were close to the magnitudes of the annual air temperature range perturbations, which generated them, reaching on average 49 % and 52 % of the perturbations for MATFS and as much as 82 % and 88 % for MATCM. I fa spanned −3309±667 to −2041±492 • C d and −3270±523 to −1873±429 • C d at the Brno-Černovice and Nebanice site, 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 290 average by as low as ∼0.6 % and ∼0.9 % of the magnitudes of the annual air temperature range perturbations.

Model performance and limitations
Generally, the model validation showed its relatively high success rate for most air temperature characteristics if appropriate driving parameters are available. It is remarkable, though, given that active-layer thickness commonly exhibits rather moder-295 ate 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 PERICLIMv1.0 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 its rather high accuracy. It also needs to be stressed that the model was successful regardless the ground physical properties used 300 have certainly undergone at least slight changes since sampling and varied over the validation period. More accurate and less scattered outputs on James Ross Island than in Alaska were largely due to differences in the distribution of the ground physical properties within the active layer, which is rather homogeneous at the James Ross Island sites (Hrbáček et al., 2017a;Hrbáček and Uxa, 2020), but has a two-layer structure (Appendix D) with a thick surface layer of peat at the Alaskan locations (Zhang, 1993). Besides, the model parameterizes the air temperature behaviour with a sine wave, which simplifies its actual evolution 305 over a year and completely ignores sub-annual variations.
Overall, however, the model underestimated most air temperature characteristics to various extents (Fig. 3). Firstly, it 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 as well as the active-layer temperature at the onset of thawing (Romanovsky and Osterkamp, 1997;Kurylyk and Hayashi, 2016 (Kurylyk and Hayashi, 2016), these require additional inputs, such as frozen thermal conductivity, thawed and frozen volumetric heat capacity, or active-layer temperature at the start of its thawing. As such, the corrections are frequently difficult to implement even in many present-day situations and definitely are much less viable for palaeo-applications. Moreover, their inverse solution would not be straightforward and would probably demand iterative techniques. Similarly, the Johansen's (1977) thermal-conductivity model is advantageous in that it requires fewer inputs as 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 another thermal-conductivity schemes.
Besides, it should be highlighted that for the model validation cases the active-layer thickness was trimmed by the depth 320 of the shallowest ground temperature sensor, which was used to determine the ground-surface thawing n-factor, to ensure consistency of the calculations (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 the effects described in the previous paragraph. Since the published ground-surface thawing n-factors used (Juliussen and Humlum, 2007;Lewkowicz et 325 al., 2012;Way and Lewkowicz, 2018) built on various depths of ground temperature sensors, we did not adjust the activelayer thickness for the palaeo-applications, and the modelled palaeo-air temperature characteristics could thus have a slightly increased accuracy as well.
By contrast, it disagrees with slightly milder MAAT reductions of at least −7 • C and −13 to −6 • C derived from groundwater data (Corcho Alvarado et al., 2011) and borehole temperature logs (Šafanda and Rajver, 2001) (Hijmans et al., 2005). The GCMs suggest that MAAT was between −4.5 • C and −0.4 • C at the study sites, which corresponds to its reduction between −12.6 • C and −9 • C. 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).

345
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 Central European lowlands (Huijzer and Vandenberghe, 1998;Marks et al., 2016). By contrast, it deviates greatly from the GCMs outputs being as high as 11.1 to 18 • C. The modelled MATCM ranged widely from −23.2±3.2 to −13.7±2.5 • C, which, however, also pertains to other proxy records that show somewhat lower values of −27 to −16.5 • C (Hui-but its average is generally well consistent with that modelled by this study. Unfortunately, other modelled palaeo-air temperature characteristics cannot be validated directly because no reliable proxy records or GCMs outputs are available for them.
However, as they are closely related to MAAT, MATWM, or MATCM, we believe that those attributes have similar plausibility.
Lastly, the modelled MAAT range between −7.0±1.9 • C and −3.2±1.5 • C is relatively well consistent with the MAAT threshold of <−4 • C commonly suggested for the formation of cryoturbation structures within fine-grained substrates (Van-355 denberghe, 2013). At the study sites, however, the cryoturbations consist of coarse-grained materials, for which the MAAT threshold as low as <−8 • C has been proposed (Vandenberghe, 2013). Although the distinction between fine-and coarsegrained substrates may not be clear, it definitely raises questions about the validity of the previously suggested MAAT thresholds for cryoturbation structures that have been widely utilized in reconstructions of past temperatures, and calls for their thorough revision.

Sensitivity analysis
Global sensitivity analysis using multiple regression and resulting standardized regression coefficients (SRCs), which indicate how many standard deviations an output parameter changes in response to one standard deviation change in an input parameter, with all other inputs being constant, 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 site (Fig. 7). The 365 palaeo-active-layer thickness importantly showed the highest 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 370 temperature range negligibly affected the thawing-season air temperature attributes (Fig. 7), and these are thus assumed to be the most plausible. Ground-surface and subsurface driving parameters, such as volumetric ground moisture content, ground dry bulk density, and ground-surface thawing n-factor, had considerably lower, albeit stable, effects on most modelled palaeo-air temperature characteristics. Ground quartz content was the weakest of the driving parameters as it was responsible only for a minor variability in the model outputs (Fig. 7).

Applications to other periglacial features
Besides cryoturbations, the model is thought to be utilized for deriving palaeo-air temperature characteristics on the basis of any other 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, mountain-top detritus, active-layer detachment slides, upfrozen clasts, indurated horizons, or frost-weathering microstructures (Ballantyne and Harris, 1994;Matsuoka, 2011;Ballan-380 tyne, 2018). However, it can also be easily adapted for seasonal-frost features.  Table 1 for abbrevations.
Since most periglacial features develop on at least decadal or centennial timescales (e.g., Karte, 1983;Matsuoka, 2001;Ballantyne, 2018), inherently involving natural climate as well as active-layer thickness variations, we hypothesize that their vertical extent probably also includes the bottom transient layer where the boundary between the active layer and permafrost has fluctuated when periglacial features formed (cf. Shur et al., 2005). This implies that the palaeo-active-layer thickness may 385 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, such as small patternedground features, small cryoturbations, some solifluction structures, up-frozen clasts, indurated horizons, or frost-weathering microstructures, may be produced by seasonal frost alone, while some look-alike features may even have a non-periglacial 390 origin (Ballantyne and Harris, 1994;Ballantyne, 2018). Nonetheless, even if identified correctly, problems may arise in places where ground-surface level has changed over time due to sedimentation or erosion. Besides, high uncertainty can probably be expected especially for periglacial features composed of pebbly to bouldery materials, such as blockfields or mountain-top detritus, in which the vertical variability of ground physical properties is extremely high (Ballantyne, 1998) that is difficult to describe by single-value parameters. On slopes, these materials may also provoke non-conductive heat-transfer processes, 395 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, the model should be ideally utilized for co-occurring periglacial features that allow for a more robust palaeo-air temperature estimates. Surely, it can capture only a short 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 exist 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 also 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., in print).

415
The model success rate is definitely promising and suggests that it could become a useful tool for reconstructing Quaternary palaeo-environments across vast areas of mid-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 and, 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 devel-420 opments of more sophisticated modelling tools that will further increase the exploitability and reliability of relict periglacial features as indicators of palaeo-climates.
Code and data availability. The latest version of PERICLIMv1.0 is available as R package from https://github.com/tomasuxa/PERICLIMv1. 0 under the GPLv3 license. The exact version of the model used to produce this paper is archived at https://doi.org/10.5281/zenodo.4057202.
The validation datasets from James Ross Island are available upon request from FH (hrbacekfilip@gmail.com), whereas those from Alaskan 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's (1977) thermal-conductivity model. This empirical transfer scheme requires less parameterizations than its derivatives, but still provides reasonable and consistent estimates of thermal conductivity for a wide range of substrates 430 having the 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 dry k dry and saturated k sat thermal conductivity (W m −1 K −1 ) of a material as follows: 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.

435
The dry ground thermal conductivity is defined by the following semi-empirical relationship (Johansen, 1977): where the constant of 2 700 kg m −3 represents the typical density of solid ground particles that is, for consistency, used throughout the thermal-conductivity scheme.
The saturated ground thermal conductivity is calculated by a weighted geometric mean based on the thermal conductivities 440 of individual ground constituents and their respective volume fractions (Johansen, 1977): where k s and k w is the thermal conductivity of solid ground particles and water, respectively (W m −1 K −1 ), and n is the ground porosity (−), which is expressed as a function of the dry bulk density and the typical density of solids: The thermal conductivity of water is set at 0.57 W m −1 K −1 , while that of solids is computed as: where k q and k o is the thermal conductivity of quartz and other minerals, respectively (W m −1 K −1 ), and q is the quartz fraction of the total content of solids (−). Quartz is assigned the thermal conductivity of 7.7 W m −1 K −1 , while for other minerals it is as follows (Johansen, 1977): Note that materials having more than 5 % of clay should be considered as 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, 455 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 the degree of saturation S (−) larger than 10 % and 5 %, respectively, is expressed as (Johansen, 1977): 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).
Appendix B: Solution using MATWM to define the range of annual air temperature oscillations 465 Unconventionally, the range of annual air temperature oscillations can also be defined using MATWM (cf. Williams, 1975)  calculated numerically using the bisection root-finding algorithm, but is searched for −∞ < MAAT ≤ 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 to estimate MATWM. However, it should be employed cautiously because 470 it is highly sensitive to variations of I ta and MATWM itself (Fig. 2), which can result in major errors if the input parameters are defined inaccurately. Unfortunately, the deviations are expected to be higher at lower both I ta and MATWM (Fig. 2) (Table C1). The stations measured air and ground temperatures with thermistor sensors installed in solar radiation shields 1.5 or 2 m above ground surface and at six to fifteen depth levels ranging from or near from the ground Table C1. Driving parameters used to model air temperature characteristics at the James Ross Island (upper section) and Alaskan Arctic (lower section) validation sites. Note that the parameters described by means and standard deviations (ξ, nt, Aa) had normal distributions, while others were assumed to be uniformly distributed (q) or were treated as constants (φ, ρ). 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).

485
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 the maximum annual ground temperature > 0 • C 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. Thawing and freezing seasons were defined by a continued prevalence of positive and negative mean daily temperatures, respectively, at the shallowest ground temperature 490 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, I ta , L t and MATFS, I fa , L f . 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 the active-layer formation than a fixed calendar period (Hrbáček 495 and Uxa, 2020). Annual air temperature range was defined by an annual spread of a 31-day simple central moving average of mean daily air temperatures, with its extremes being considered to substitute MATWM and MATCM. Finally, ground-surface thawing n-factor was derived as a ratio of the thawing index at the shallowest ground temperature sensor and air thawing index.
Consequently, for modelling, the active-layer thickness had to be trimmed by the depth of the shallowest ground temperature sensor to ensure consistency because the model presumes that the ground-surface thawing n-factors transfer between air and 500 ground-surface temperatures (cf. 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/2014 to 2018/2019 (Hrbáček et al., 2017a;Hrbáček and Uxa, 2020), while those for the Alaskan sites were adapted from Zhang (1993) and Romanovsky and 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 a wet sieving and X-ray diffraction (Hrbáček et al., 2017a;Hrbáček and Uxa, 2020) or assessed visually via soil type (Zhang,510 1993). All the substrates were considered as 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 activelayer thickness in two-layer ground can be applied. It has been proposed in the following form (Nixon and McRoberts, 1973;Kurylyk, 2015): where Z 1 is the thickness of the top sub-layer (m), the physical parameters of which are subscripted by 1, while the bottom sub-layer is denoted by the subscripts 2. The ground-surface index can then be simply expressed from Eq. (D1): Lφ 2 ρ w 2k t2 .
As in Eq. (1 and 2), the product of φ and ρ w can be substituted by that of the gravimetric moisture content and dry bulk density 520 of the ground, but note that the fraction on the far right of Eq. (D1) and at the corresponding place of 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 (Eq. 3 to 14).
Author contributions. TU came up with an initial idea with feedbacks from MK, developed the model, evaluated it against data from James Ross Island and Alaskan Arctic, which were processed by FH and TU, respectively, and tested it for derivation of palaeo-air temperature characteristics using data sampled collectively with MK in the Czech Republic. TU draw figures and wrote the manuscript 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.