Articles | Volume 17, issue 6
Research article
10 Dec 2021
Research article |  | 10 Dec 2021

A new perspective on permafrost boundaries in France during the Last Glacial Maximum

Kim H. Stadelmaier, Patrick Ludwig, Pascal Bertran, Pierre Antoine, Xiaoxu Shi, Gerrit Lohmann, and Joaquim G. Pinto

During the Last Glacial Maximum (LGM), a very cold and dry period around 26.5–19 kyr BP, permafrost was widespread across Europe. In this work, we explore the possible benefit of using regional climate model data to improve the permafrost representation in France, decipher how the atmospheric circulation affects the permafrost boundaries in the models, and test the role of ground thermal contraction cracking in wedge development during the LGM. With these aims, criteria for possible thermal contraction cracking of the ground are applied to climate model data for the first time. Our results show that the permafrost extent and ground cracking regions deviate from proxy evidence when the simulated large-scale circulation in both global and regional climate models favours prevailing westerly winds. A colder and, with regard to proxy data, more realistic version of the LGM climate is achieved given more frequent easterly winds conditions. Given the appropriate forcing, an added value of the regional climate model simulation can be achieved in representing permafrost and ground thermal contraction cracking. Furthermore, the model data provide evidence that thermal contraction cracking occurred in Europe during the LGM in a wide latitudinal band south of the probable permafrost border, in agreement with field data analysis. This enables the reconsideration of the role of sand-wedge casts to identify past permafrost regions.

1 Introduction

Permafrost is an important component of the climate system and is particularly sensitive to variations in climate. Permafrost is defined as ground – including soil, rock, ice, and organic material – that remains at or below 0 C for at least 2 consecutive years (e.g. van Everdingen2005). In recent decades, thawing of permafrost soils has affected many high-latitude regions, and thawing is most likely to accelerate in the near- and long-term future (IPCC2013, 2019; Harris et al.2009). While enhanced greenhouse gas forcing leads to warming temperatures and, thus, to permafrost thawing, the thawing itself leads to the release of greenhouse gases that were previously bound within the frozen soils. Therefore, the greenhouse effect is enhanced and leads to further warming of the climate in a positive feedback (e.g. IPCC2019; Liu and Jiang2016a; Schuur et al.2015).

Current climate model simulations project a large range of uncertainties regarding the decrease in permafrost areas (e.g. IPCC2019; Schuur et al.2015). The models are calibrated for present-day conditions, under which they are well tested. However, the responses of the models to the same forcing vary by several orders of magnitude (e.g. Braconnot et al.2012; Cleator et al.2020; IPCC2013, 2019). It is therefore necessary to evaluate the climate models under a wider range of climate conditions. This can be achieved by simulating past climates and comparing the results with proxy evidence from Quaternary sequences (Braconnot et al.2012; Harrison et al.2015; Smerdon2012).

During the Last Glacial Maximum (LGM; Clark et al.2009; Mix et al.2001), corresponding to around 26.5–19 kyr BP, huge ice sheets covered large parts of the Northern Hemisphere, modifying the surface albedo and orography (Hughes et al.2015; Ullman et al.2014), and enhanced sea ice cover modified heat fluxes between the ocean and atmosphere (Flückiger et al.2008). During the coldest phase of the LGM, the sea level was about 130 m lower than today (Lambeck et al.2014) and the greenhouse gas concentrations were at a historical minimum with values less than half of present-day concentrations (Clark et al.2009; Monnin et al.2001). Lower greenhouse gas concentrations favoured the growth of C4 over C3 plants (Prentice and Harrison2009), although only C3 plants have actually been identified in European loess (Hatté et al.1998, 2001). Globally, this hampered the development of trees (Woillez et al.2011), resulting in less-productive terrestrial ecosystems and more open vegetation (Bartlein et al.2011). Ultimately, this induced easily erodible soils, whose contribution to the dust cycle increased (Prospero et al.2002; Ray and Adams2001). These boundary conditions and forcing led to a substantially different climate than today. In general, the LGM was a colder, drier, and windier period in Earths' history compared with the recent climate (e.g. Annan and Hargreaves2013; Bartlein et al.2011; Löfverström et al.2014). The global and annual mean surface air temperatures were about 4 C colder than today, with differences reaching up to 14 C close to the LGM ice sheets in areas such as central Europe (e.g. Annan and Hargreaves2013; Bartlein et al.2011; Clark et al.2009; Pfahl et al.2015; Ludwig et al.2017). The atmospheric circulation in the North Atlantic region varied considerably from the current conditions, mainly due to the direct influence of the altered topography by ice sheets (Justino and Peltier2005; Merz et al.2015). A planetary large-scale atmospheric wave with an amplitude much larger than today was induced, with a deep trough downstream of the Laurentide ice sheet. This led to a generally more zonal orientation of the North Atlantic jet stream (Löfverström et al.2014). Additionally, the jet was enhanced and its position was shifted southward (e.g. Li and Battisti2008; Merz et al.2015; Pausata et al.2011). The storm track during the LGM evolved accordingly (e.g. Löfverström et al.2014; Ludwig et al.2016; Raible et al.2021), and extreme cyclones were more intense and characterised by less precipitation (Pinto and Ludwig2020). Thus, cyclones were able to trigger more frequent dust storms during the LGM (Antoine et al.2009; Pinto and Ludwig2020; Sima et al.2009). Besides these dust storms, easterly winds induced by an anticyclone over the Fennoscandian ice sheet (FIS) were another important factor for the deposition of loess in central and western Europe (Krauß et al.2016; Schaffernicht et al.2020; Stevens et al.2020) as well as westerly to north-westerly winds (e.g. Renssen et al.2007; Schwan1986, 1988).

At the same time, adjacent areas south of the FIS were widely affected by permafrost (Kitover et al.2013; Levavasseur et al.2011; Saito et al.2013; Vandenberghe et al.2014; Washburn1979). The past permafrost distribution is usually inferred from the occurrence of a variety of fossil periglacial features, among which ice-wedge pseudomorphs are the most reliable and widespread (e.g. Bertran et al.2014; Huijzer and Isarin1997; Péwé1966; Vandenberghe1983; Vandenberghe et al.2014). Ice wedges develop within perennially frozen ground, when the temperature drops quickly and the ground experiences thermal contraction cracking. Annual frost cracks that reach downward into the permafrost are a few millimetres wide. They get filled with snowmelt that freezes into ice veins. Repeated cracking over years at the same location adds ice veins that constitute ice wedges (e.g. Harry and Gozdzik1988; Murton2013). Ice-wedge pseudomorphs observed from the LGM in Europe were formed when the ice melted and the cavities were filled by collapsing soil materials. Today, ice wedges are mostly active in continuous permafrost environments (Fortier and Allard2005; Kokelj et al.2014; Matsuoka et al.2018; Péwé1966). Open cracks may also be filled with wind-blown sand, which gives rise to sand wedges, or by both ice and sand, which gives rise to composite wedges. Active sand wedges are currently primarily found in areas characterised by continuous permafrost and limited snow and vegetation cover (i.e. the polar deserts), and with local sources of aeolian sediments, such as in Antarctica (Bockheim et al.2009; Levy et al.2008; Murton et al.2000; Péwé1959). Ground cracking is often restricted to the active layer (i.e. the surface layer subjected to seasonal freezing and thawing) in the areas underlain by “warm” permafrost (i.e. at a temperature close to 0 C) and south of the permafrost border. Thin cracks develop and are referred to as seasonal frost cracks. However, Wolfe et al. (2018) showed that large shallow sand wedges can also develop in Canada in areas with deep seasonal ground freezing (i.e. without perennially frozen ground) in mineral soils close to dune fields, which provide abundant sand to fill the cracks.

Thermal contraction cracking of the ground is the causal factor that leads to ice (or sand) wedge growth. Ecological factors such as type of vegetation cover and thick snow cover often limit thermal contraction cracking, as they may prevent the cooling of the ground. This is the case in current densely vegetated areas that insulate the ground and trap snow (e.g. shrub tundra and taiga; Kokelj et al.2014; Mackay and Burn2002). Conversely, cracking can occur at low frequency in mid-latitude, cool temperate regions in grounds devoid of tall vegetation and snow, particularly in roads and airport runways (Barosh2000; Okkonen et al.2020; Washburn1963).

Many attempts at reconstructing the past permafrost distribution in Europe using field proxies have been performed during the last decades. Based on the assumption that both active ice wedges and sand wedges are associated with continuous permafrost and possibly with widespread discontinuous permafrost (Burn1990; Romanovskij1973), some of the earliest reconstitutions, as reported by Vandenberghe et al. (2014), proposed that Europe was affected by permafrost as far south as 43.5 N. However, a detailed analysis of periglacial features in France by Andrieux et al. (2016b, 2018) demonstrated that typical ice-wedge pseudomorphs are exclusively found north of 47.5 N, whereas sand-wedge casts occur at lower latitude at the periphery of aeolian sand sheets. A correlation between wedge depth and latitude has also been highlighted, which strongly suggests that the southernmost shallow sand wedges developed in regions where perennial ice could not form, i.e. without permafrost or with sporadic permafrost. A similar pattern has also been highlighted in China by Vandenberghe et al. (2019). The sand wedges reach up to 1 m wide in south-west France near 45 N in the periphery of cover sands. Optically stimulated luminescence dating of the sand fill by Andrieux et al. (2018) demonstrated that these large epigenetic sand wedges resulted from repeated periods of growth throughout the Last Glacial.

Multiple attempts have also been performed to infer the LGM permafrost occurrence from climate model data. Liu and Jiang (2016b) considered both direct and indirect methods. The simplest indirect method is based on the modelled mean annual air temperature (MAAT). Threshold values for permafrost occurrence were adapted according to ground texture (Vandenberghe et al.2012). However, this method only provides a rough estimate of permafrost extension, as a variety of other factors are known to impact ground temperatures, including water content, vegetation, and snow cover. Particularly, the insulating effects of snow and vegetation cover may be responsible for an offset of up to 6 C between the MAAT and the mean annual ground surface temperature (MAGST). On the other hand, variations in ground thermal conductivity (depending on texture and water content) may result in an offset of 2 C between the MAGST and the temperature at the top of permafrost (TTOP) (e.g. Smith and Riseborough2002; Throop et al.2012). A refined indirect method to derive permafrost occurrence from climate model data is the use of the surface frost index (SFI, Nelson and Outcalt1987), which corresponds to the ratio between frost and thaw penetration depths and takes the effects of snow in account. The SFI has been used in several studies, with only minor changes to the original method. For example, monthly model output was used instead of summing up daily air temperatures (e.g. Frauenfeld et al.2007; Liu and Jiang2016b). Slater and Lawrence (2013) weighted the snow depth for each month to consider snow accumulation effects, whereas Stendel and Christensen (2002) replaced the surface air temperature with the temperature of their deepest simulated ground layer (5.7 m deep) to investigate permafrost degradation due to current global warming. The latter authors pointed out the advantage of taking simulated ground temperatures, where insulation effects of snow and vegetation cover are explicitly taken into account by the models, and rendered empirical approaches redundant. For the direct method, the modelled ground temperatures below 0 C are used to diagnose permafrost. The studies differ slightly with respect to the depth of the considered ground temperatures (e.g. Liu and Jiang2016a, b; Saito et al.2013; Slater and Lawrence2013).

Studies investigating the permafrost limits during the LGM using global climate simulations have so far failed to appropriately reproduce the permafrost extent as reconstructed from field proxies (e.g. Andrieux et al.2016b; Levavasseur et al.2011; Ludwig et al.2017). However, there is evidence for improvements when using the data from regional climate simulations (e.g. Ludwig et al.2017, 2019). The aim of this study is (1) to explore the possible benefit of using regional climate model data to improve the permafrost representation over France, (2) to decipher how the atmospheric circulation affect the permafrost boundaries in the models and finally, (3) to test the role of ground thermal contraction cracking in wedge development during the LGM.

In Sect. 2, we introduce the adaptions made to the regional climate model to be compliant with LGM boundary conditions and describe the global simulations that provide the initial and boundary conditions. Further, we give an overview of the different methods used to derive the LGM permafrost distribution in France. In Sect. 3, we describe the general characteristics and differences of the LGM climate based on the global and regional climate model data and present the permafrost and ground cracking distribution based on regional climate model data. Finally, we discuss and summarise the results in Sect. 4.

2 Data and methods

In this study, LGM simulations of two global climate models, namely MPI-ESM-P (MPI – Max Planck Institute; Jungclaus et al.2013; Stevens et al.2013) and AWI-ESM (AWI – Alfred Wegener Institute; Sidorenko et al.2015; Lohmann et al.2020), are dynamically downscaled with the Weather Research and Forecasting model (WRF; Skamarock et al.2008). Both global models share the same atmospheric component ECHAM6 but different modules for the ocean. The MPIOM (Marsland et al.2003) is coupled within the MPI-ESM-P, forming the well-established global climate model that took part in several Coupled Model Intercomparison Project (CMIP) phases. In the AWI-ESM, the FESOM ocean model (Wang et al.2014) featuring an unstructured mesh as well as a multi-resolution approach is used with a relatively high resolution of less than 30 km north of 50 N. The atmospheric grid applied in the MPI and AWI experiments is T63 (roughly 1.9 spatially) with 47 unevenly distributed vertical levels. The simulations follow either the Paleoclimate Modelling Intercomparison Project Phase 3 (PMIP3) protocol (MPI; Braconnot et al.2012;, last access: 6 December 2021) or the PMIP4 protocol (AWI; Kageyama et al.2017), where the boundary conditions (solar constant, orbital parameters, greenhouse gases) are set according to the best estimate of the LGM boundary conditions. The AWI-ESM has been used in the recent CMIP6/PMIP4 intercomparisons (Brierley et al.2020; Keeble et al.2021; Kageyama et al.2021) and was applied for the LGM (Lohmann et al.2020). The ice sheet provided for PMIP3/CMIP5 LGM experiments is a blended product obtained by averaging three different ice sheet reconstructions: ICE-6G v2.0 (Peltier et al.2015), MOCA (Tarasov and Peltier2003), and ANU (Lambeck et al.2002). In contrast, the LGM topography in the AWI experiment is configured based on the ICE6G reconstruction (Peltier et al.2015). For the recent climate, the pre-industrial period (PI), corresponding to roughly 1850, is used as a reference. The simulations again follow the PMIP3 (MPI; Taylor et al.2012) or PMIP4 protocol (AWI; Eyring et al.2016).

To account for model uncertainties, outputs from these global LGM simulations are used to drive the regional WRF simulations. The atmospheric boundary conditions are updated every 6 h, and sea surface temperature (SST) and sea ice cover are updated daily. Apart from the different forcing, the set up of the two regional simulations is identical. The coastlines, ice sheet extent, trace gas conditions, and orbital parameters are adapted to LGM values according to the PMIP3 protocol (Ludwig et al.2017). Modifications to the Alpine ice sheet are implemented according to Seguinot et al. (2018). Land use and vegetation cover is taken from the CLIMAP data set (CLIMAP Project Members1984). An overview of the parameterisation schemes used in the WRF simulations is given in Table 1. Most important for the representation of the ground characteristics is the parameterisation of the land surface, for which we used the unified Noah land surface model (Tewari et al.2004). Based on 19 different soil types, various ground parameters (e.g. ground thermal conductivity) are set and used for the calculations of ground temperatures and moisture for each grid point. More details can be found in studies such as Chen and Dudhia (2001) and Niu et al. (2011) as well as references therein. The first model domain covers large parts of Europe with a horizontal resolution of 50 km (see Fig. 1) and 35 vertical layers up to 150 hPa. The integration time step is 240 s. The second, nested domain covers southern parts of the FIS, the Alps, and France, where the latter represents the target region to assess the LGM permafrost limits in this study. Here, the horizontal resolution is 12.5 km and the integration time step is 48 s. The soil is separated into four layers, with representative depths of 5, 25, 70, and 150 cm. A total of 32 years are simulated for each global forcing simulation. The first 2 years are used as a spin-up phase and are excluded from further analysis. Thus, it is ensured that the atmosphere and soil properties and processes are in equilibrium.

Hong et al. (2004)Iacono et al. (2008)Jiménez et al. (2012)Tewari et al. (2004)Hong et al. (2006)Kain (2004)

Table 1Physical parameterisation schemes used for the regional simulations.

Download Print Version | Download XLSX

Figure 1WRF model domain orography with the LGM coastline (black line) and LGM ice sheet (pink line). (a) Domain 1 with 50 km grid spacing; (b) domain 2 with 12.5 km grid spacing. Locations A and B on the map refer to the time series of soil temperatures and soil temperature gradients shown in Fig. 7 and Fig. S6 in the Supplement.

The permafrost distribution is derived from climate model data using the three different methods described in Sect. 1. For MAAT, the 2 m air temperature is considered. Threshold values were derived from data compiled from studies in current Arctic regions, where continuous permafrost is inferred for MAATs <-8±2C, whereas discontinuous permafrost requires MAATs <-4±2C (e.g. Smith and Riseborough2002; Vandenberghe et al.2012). The surface frost index (SFI) is based on the annual freezing and thawing degree-days (DDF and DDT respectively) which refer to the sum of daily air temperatures below or above 0 C respectively:


An SFI between 0.5 and 0.6 indicates sporadic permafrost, between 0.6 and 0.67 indicates discontinuous permafrost, and above 0.67 indicates continuous permafrost (e.g. Nelson and Outcalt1987; Stendel and Christensen2002). For this indirect method, we use ground temperatures of the third layer at 78 and 70 cm for the global and regional simulations respectively. With the direct method, permafrost is inferred when ground temperatures are at or below 0 C.

Beyond the permafrost indices, ground cracking is assumed to be possible when two conditions derived from fieldwork by Matsuoka et al. (2018) are fulfilled simultaneously: a daily mean soil temperature below 5 C at a depth of 1 m and a temperature gradient in the upper metre of the ground below 7 Cm-1. These minimum values might represent shallow cracking within the active layer or seasonally frozen layer and can be compared against the sand-wedge distribution. Conditions for intensive and deep thermal contraction cracking (T100=-10C and GAL=-10Cm-1) are tested in regard to the ice-wedge pseudomorph distribution in France. Due to higher ice content and higher organic carbon content of the ground, these values do not necessarily correspond exactly to those of France during the Pleistocene. We use the third soil layer again, with depths of 78 cm in the global simulations and of 70 cm in the regional simulations.

To evaluate the model simulations, the distribution of ice-wedge pseudomorphs and sand wedges after Andrieux et al. (2016b) and Isarin et al. (1998) are considered.

3 Results

3.1 Global boundary conditions

In this section, we present the large-scale characteristics of the LGM climate derived from global climate model data that is used for dynamical downscaling in comparison with the respective PI simulations. It is important to investigate the climatic mean state and possible biases of the global projections in order to be able to interpret the regional simulations accurately.

Both global models simulate colder annual mean SSTs under LGM than under PI conditions (see Fig. 2a and b). For the MPI model, a limited area with enhanced SSTs is simulated over the North Atlantic. This does not match with proxy data (MARGO Project Members2009) and is a known issue for this and other PMIP3 models (e.g. Wang et al.2013; Ludwig et al.2016, 2017). The AWI simulation does not show this warm anomaly over the North Atlantic and the SSTs are generally colder. In the Arctic Ocean, the SSTs in the AWI simulation are considerably higher than in the MPI simulation. This can be explained by the sea ice cover, which is lower in the AWI LGM simulation.

Figure 2(a) Distribution of annual mean SST differences between the global MPI simulations under LGM and PI conditions. Panel (b) is the same as panel (a) but for the global AWI simulations. Panels (c) and (d) are the same as panels (a) and (b) but for the annual mean wind speed at 300 hPa.

The analysis of wind speed at 300 hPa gives insights into the jet stream structure and strength, which are dominant factors of the atmospheric large-scale circulation over the North Atlantic/European region. In agreement with Li and Battisti (2008), both models show a stronger jet under LGM conditions compared with the simulations under PI conditions (see Fig. 2c and d). This is particularly the case over the North Atlantic, south-eastward of the Laurentide ice sheet, where the annual mean wind speed is up to 14 m s−1 higher for the LGM. On the other hand, the wind speed on both the southern and northern flanks of the jet stream is actually 2–4 m s−1 weaker during the LGM, indicating a more constrained large-scale flow. Even though the wind anomalies are quite similar for both global climate models (GCMs), the actual structure is dissimilar: while the jet is less constrained and deflected to the north for the MPI simulations, reaching Europe at the latitude of Ireland, the jet stream in the AWI GCM reaches Europe at the latitude of the Iberian Peninsula and France and extends farther into the continent. In general, the simulated winds speeds at 300 hPa in the AWI model are weaker compared with the MPI model (not shown). The zonal structure of the wind speed anomalies identified for the AWI simulations is more similar to the ensemble mean of CMIP5 models (e.g. Ludwig et al.2016) than the MPI anomaly pattern.

3.2 Climate of the regional simulations

Based on the GCM simulations, we obtain two different variants of the regional LGM climate in western Europe. The results are shown primarily for the larger domain of the regional simulations, as the climate of the high-resolution simulations yields a similar structure.

The annual mean 2 m air temperature is considerably lower in the WRF-MPI than in the WRF-AWI simulation (see Fig. 3). The biggest differences are identified near the ice sheet margin – almost 10 C in the respective annual means. The sign and pattern of the differences are visible in all seasons, but winter air temperatures clearly diverge most. For the summer, air temperatures in both models are more similar to each other. These differences can be partly attributed to the snow cover: except for summer, almost the entire region is covered by snow in WRF-MPI, even though a snow height of several metres is only reached over the FIS and the Alpine region (see Fig. S3 in the Supplement). WRF-AWI shows markedly higher snow accumulation over the ice sheets with differences of more than 20 m compared with the WRF-MPI simulation but generally shows less snow cover in southern and central Europe. Differences amount to 20 % less snow cover in WRF-AWI in the respective annual means and to 40 % in both spring and winter. In summer, only the ice sheets are snow covered in both simulations; thus, the differences are negligible. Nevertheless, more precipitation is simulated over Europe in the WRF-AWI simulation (see Fig. 4). High precipitation amounts are either orographically induced, as for precipitation over the Alps and over the FIS, or they are associated with the moisture availability of the North Atlantic.

Figure 3Distribution of 2 m air temperature in annual (a–c) and seasonal winter (d–f) and summer (g–i) means as simulated with the regional WRF model with MPI forcing (a, d, and g) and with AWI forcing (b, e, and h) as well as their differences (c, f, and i). The black line shows the LGM coastline, and the pink line denotes the LGM ice sheet.

Figure 4Distribution of total precipitation in annual (a–c) and seasonal winter (d–f) and summer (g–i) means as simulated with the regional WRF model with MPI forcing (a, d, and g) and with AWI forcing (b, e, and h) as well as their differences (c, f, and i). The black line shows the LGM coastline, and the pink line denotes the LGM ice sheet.

Figure 5Distribution of 10 m wind speed in annual (a–c) and seasonal winter (d–f) and summer (g–i) means as simulated with the regional WRF model with MPI forcing (a, d, and g) and with AWI forcing (b, e, and h) as well as their differences (c, f, and i). The black line shows the LGM coastline, and the pink line denotes the LGM ice sheet.

The absolute annual mean wind field and the associated differences are depicted in Fig. 5. Both simulations show strongest winds south of the FIS in the respective annual and winter means, although with a notably enhanced pattern in WRF-MPI, where this also holds for each season. These winds are easterlies/north-easterlies. In contrast, westerly winds from the North Atlantic are stronger in WRF-AWI and, thus, transport heat and moisture towards Europe. During winter, the westerly winds are directed towards the centre of the domain in WRF-AWI, whereas the winds have a more south-western component in WRF-MPI and are directed towards the outside of the domain. In summer, both the WRF-MPI and WRF-AWI simulations are characterised by westerly winds from the North Atlantic. Again, winds from the FIS are blowing south- and south-eastwards, but the summer wind speeds are consistently weaker than in winter for both simulations. These wind fields are induced by the large-scale circulation in the global forcing simulations. In fact, the northerly and easterly components predominantly occur in the MPI simulation (Ludwig et al.2016), whereas southerly and westerly components occur more often in the AWI simulation. This is in accordance with the jet structure in both global simulations. As the influence of the ice sheet is higher in the (global and regional) MPI simulations, this is consistent with a partially drier and generally colder climate in western Europe during the LGM.

3.3 Permafrost and ground cracking distribution

The permafrost distribution of the global and regional simulations based on the SFI is depicted in Fig. 6. The permafrost extent based on the AWI-ESM and WRF-AWI simulations does not reach farther south than the ice sheet, apart from the Alps in WRF-AWI. A modest increase in the permafrost area is simulated by the global MPI simulation. Here, continuous permafrost is still limited to the ice sheet, but sporadic permafrost is slightly more widespread. The WRF-MPI simulation shows a larger permafrost extent. In eastern Europe, the distribution of ice-wedge pseudomorphs (Isarin et al.1998) strictly overlaps that of modelled continuous permafrost in the selected layer with a depth of 70 cm. In western Europe, field evidence for permafrost exceeds the modelled sporadic permafrost to the south. The conditions for discontinuous and sporadic permafrost are rarely fulfilled in all simulations.

Figure 6Permafrost distribution based on the surface frost index (SFI) at approximately 70 cm depth as simulated by the global MPI (a) and the global AWI (b) simulations and their respective regional counterpart (c and d). Ice-wedge pseudomorphs, composite wedges, and sand wedges from Andrieux et al. (2016) and Isarin et al. (1998) are denoted by the purple, grey, and red triangles respectively. The black line shows the LGM coastline, and the pink line denotes the LGM ice sheet.

The results of the direct method (see Fig. S4 in the Supplement) using long-term mean annual soil temperatures agree with the permafrost extent based on the SFI. However, the different types of permafrost cannot be distinguished by this method, leading to a permafrost line that corresponds to that of the sporadic permafrost based on the SFI.

Permafrost estimations based on MAAT are limited to the permanent ice areas during the LGM in all four simulations (see Fig. S5 in the Supplement). Despite the different regional climates, the reconstructed permafrost boundaries in this study closely resemble each other for MAAT. The regional climate model simulations show some additional permafrost areas, which are related to higher orography, especially in the Alps, and, in WRF-MPI, also in the Pyrenees and the Massif Central (see Fig. 1 and Fig. S2 in the Supplement). These mountainous areas are not adequately resolved in the global forcing simulations because of the coarse horizontal grid spacing.

Conditions for thermal contraction cracking after Matsuoka et al. (2018) have been tested based on the global and regional climate model data. Examples of how the soil temperature and the gradient develop over 2 consecutive years in France (locations A and B in Fig. 1) are shown in Fig. 7. Time series of the entire simulation periods for these locations can be found in Fig. S6. The two minimum criteria, (a) ground temperature at 1 m below 5 C and (b) temperature gradient in the upper metre of ground greater than 7Cm-1, are fulfilled when both curves reach below the depicted reference line. In the WRF-MPI simulation (Fig. 7a, c), this is the case several times in both years and locations, but it is not the case in the WRF-AWI simulation.

Figure 7Time series of the daily mean soil temperature and soil temperature gradient at locations A (model grid point 49.8 N, 2.49 E) and B (model grid point 46.1 N, 0.33 E), as denoted in Fig. 1, simulated in WRF-MPI (a and c) and in WRF-AWI (b and d) for 2 consecutive years. Blue lines show the development of the soil temperatures in layer 3 (70 cm). When the temperatures fall below 5 C, the first condition for thermal contraction cracking after Matsuoka et al. (2018) is fulfilled, marked with blue shading and the reference line. The soil temperature gradient between the first layer (5 cm depth) and the third layer (grey shading) is only depicted when condition two after Matsuoka et al. (2018) is fulfilled, with a gradient below 7 Cm-1. Red lines indicate the coincidence of the two conditions.


Figure 8Heat maps of the mean number of days per year when the minimum conditions for thermal contraction cracking after Matsuoka et al. (2018) are fulfilled for each grid box in the global MPI (a) and AWI simulations (b), for the first domain of the regional WRF-MPI (c) and WRF-AWI (d) simulations, and for the second domain in WRF-MPI (e) and in WRF-AWI (f). Ice-wedge pseudomorphs and sand wedges from Andrieux et al. (2016) are highlighted using cyan and red triangles respectively, only when located in France. The black line is the LGM coastline, and the grey line denotes the LGM ice sheet.


For each grid cell, the number of days per year when the thermal contraction cracking criteria (Matsuoka et al.2018) are fulfilled is translated into heat maps for each simulation. The results of the minimum conditions for (shallow) cracking are shown in Fig. 8, and the results of the conditions for intensive and deep cracking are shown in Fig. S7 in the Supplement. While the permafrost area is much smaller in the global models than their respective regional counterpart, the opposite is the case for thermal contraction cracking areas. The global AWI simulation almost meets the boundaries of sand-wedge occurrence. According to the global MPI simulation, thermal contraction cracking would have been possible as far south as the Iberian Peninsula, where no field evidence for it has been found so far. This can be associated with the lower resolution of the global simulations. Here, the Pyrenees are not resolved adequately in the model and do not act as a natural barrier for cold air arriving from the North, which can, thus, reach further south in the GCMs (see Fig. 1 and Fig. S2 in the Supplement). As for the permafrost distribution, the possible thermal contraction cracking occurrence is also poorly represented in WRF-AWI and is not able to explain the occurrence of wedges in middle and southern France. By contrast, the WRF-MPI simulation agrees well with proxy evidence. Apart from two sand wedges in the lower Rhône valley, the conditions for thermal contraction cracking are found in the simulation in the area where the features are found. This spatial coherence is further improved in the high-resolution simulation (see Fig. 8e), which can be primarily attributed to a higher resolved orography (see Fig. 1b). Moreover, the conditions for deep ground cracking are represented best in the regional WRF-MPI simulation. The heat maps show that those conditions did not occur in south-western France, which is in agreement with the field data. In this area, the sand-wedge casts do not exceed a depth of 2 m and ice-wedge pseudomorphs are not mapped at all (Andrieux et al.2016b).

4 Summary and discussion

In this study, we explore the benefit of using regional climate model data for the delimitation of the LGM permafrost distribution in comparison with field proxies in France. The main findings can be summarised as follows:

  1. The SFI is suitable to infer LGM permafrost from climate model data. The results based on the SFI are supported by the direct method, as the boundaries between permafrost occurrence and absence, as indicated by the SFI, fully match the permafrost border derived from the annual mean ground temperature. Among the models used, the SFI-based permafrost extent of the regional WRF-MPI simulation best agrees with proxy data and is clearly improved compared with its global counterpart.

  2. The thermal contraction cracking may have occurred much further south than the simulated permafrost limits, in a context of low and sparse vegetation. The southern extent of sand wedges and that of ice-wedge pseudomorphs in France as delineated by Andrieux et al. (2016b) fit well with the boundaries of LGM thermal contraction cracking derived from the regional WRF-MPI simulation based on the criteria for shallow and for deep cracking after Matsuoka et al. (2018) respectively. In contrast, the global MPI simulation does not resolve orographic features (e.g. the Pyrenees and the Rhône Valley) sufficiently, leading to a possible southward airflow transporting cold air across France to Spain, and allows ground cracking to occur at excessively low latitudes.

  3. The obtained estimates for the possible location of permafrost is consistent with the hypothesis proposed by Andrieux et al. (2016b, 2018), who suggest that sand wedges did not exclusively form in permafrost areas during the LGM but also developed within deep seasonally frozen ground. Contrary to what occurs today in large Arctic areas underlain by permafrost, where ground insulation is limited by dense vegetation (shrub tundra, taiga) and snow cover prevents ground cracking and limits the growth of ice wedges (existing ice wedges that have formed in relation to different climatic or ecological conditions do not melt but are dormant), ice-wedge growth in permafrost areas in France during the LGM was rapid because thermal conditions leading to ground cracking occurred with high frequency. Large ice wedges (which after thawing developed into recognisable pseudomorphs) would have formed in permafrost where it was cold enough in winter to crack. Simulations show that periods of winter ground temperatures below 10 C at 1 m depth could occur in the discontinuous and sporadic permafrost zone, suggesting that thermal contraction cracks were possibly not restricted to the active layer but could propagate into the permafrost in these areas leading to the development of ice wedges. The regional WRF-MPI simulations best match the proxy-based permafrost reconstruction. The agreement with the proxies is better in eastern Europe, even though the availability of field data remains scarce in that region compared with western Europe. The presence of ice-wedge pseudomorphs in northern France actually shows that permafrost must have extended at least 150 km further south than simulations suggest.

The consideration of two global models enables the quantification of uncertainties associated with the large-scale flow under LGM conditions. The global MPI and AWI simulations differ in their atmospheric flow and jet structure. In the AWI, the westerly flow dominates so that moisture and heat are transported from the North Atlantic towards Europe. This large-scale circulation is in good agreement with the multi-model mean of the CMIP5/PMIP3 and CMIP6/PMIP4 models, whereas the MPI simulation exhibits a more northward jet stream and suggests a stronger ice sheet influence through prevailing north- and north-easterly winds (Kageyama et al.2021; Ludwig et al.2016). Considering that the regional WRF-MPI simulation is largely in agreement with proxy evidence for both the permafrost and ground cracking extent, we assume that the large-scale circulation of the LGM is reflected more accurately in this simulation. For wind and air pressure, only indirect proxy evidence currently exists, e.g. the reconstruction of easterly wind directions from sediments across the European loess belt (Dietrich and Seelos2010; Krauß et al.2016; Römer et al.2016). Because of the drier conditions with less vegetation and higher wind speeds, dust events occurred frequently during the LGM. This is reflected by the thick loess deposits in western and central Europe, which form the European less belt (e.g. Lehmkuhl et al.2016). Recent studies similarly support the hypothesis that, besides individual cyclone events, easterly winds induced by a semi-permanent anticyclone over the FIS were an important component for the glacial dust cycle (e.g. Raible et al.2021; Schaffernicht et al.2020; Stevens et al.2020).

Overall, the new regional climate simulations largely reconcile the field data and enable the reconsideration of the significance of ice-wedge pseudomorphs and sand-wedge casts for understanding past climate variations. Field data still suggest a wider extension of permafrost in western Europe than shown by the simulations; however, analysing the southern extent of thermal contraction cracking completes the picture. Various factors may account for a remaining gap between proxy and model data. These factors include the following:

  1. The ground thermal conductivities used in the models may not be perfectly adequate. For fine-grained soils such as loess (in which many ice-wedge pseudomorphs have been reported), this could lead to a slightly colder ground temperature, although this effect is assumed to have been minor.

  2. Snow depth and snowpack properties (e.g. Royer et al.2021) are very sensitive factors for permafrost, and some snow processes are not considered in the models. This may explain some of the discrepancies between field data and simulations. Snow sweeping by the wind at some sites, especially on plateaus, may have led to local permafrost development. However, it should be mentioned that pseudomorphs have been described in the Last Glacial floodplains in the Paris Basin (e.g. Bertran et al.2018), i.e. in places that are favourable to snow accumulation a priori.

  3. Data from loess sections in northern France (Antoine et al.2003, 2014) and Germany (Meszner et al.2013) show that the main phases of ice-wedge development occurred between 30 and 24 ka. This period, called the Last Permafrost Maximum (LPM, Vandenberghe et al.2014), covers short and very cold events, which resulted in wider permafrost extension than during the LGM sensu stricto. However, boundary conditions for the simulations are only known accurately at 21 ka.

To conclude, the combination of the well-established permafrost index SFI and the criteria for thermal contraction cracking by Matsuoka et al. (2018), both based on regional climate model data, provides new possibilities for the estimation of the permafrost extent and the interpretation of ice and sand wedges, especially for palaeoclimate applications. In this context, the use of regional climate model simulations with a highly resolved orography is clearly beneficial (e.g. Ludwig et al.2019) and should be considered for regions other than western Europe.

Code and data availability

The source code of AWI-ESM is available from the AWI-based SVN repository (, FESOM2021). The data are available in Kageyama et al. (2021). Data used in the MPI-ESM-P PI and LGM simulations are available in Jungclaus et al. (2012a) (DOI: and Jungclaus et al. (2012b) (DOI: respectively. The WRF data will be archived at DKRZ (German Climate Computing Centre) and are available from the corresponding author upon request. PMIP3 boundary conditions can be obtained at (last access: 6 December 2021). Vegetation cover and land use data from CLIMAP (1984) can be obtained at (last access: 6 December 2021) (IRI2015). The database of ice-wedge pseudomorphs and sand-wedge casts in France is available at and discussed in Andrieux et al. (2016a). The database of European ice-wedge pseudomorphs and sand-wedge casts is available at (Huijzer et al.1998).


The supplement related to this article is available online at:

Author contributions

PL, PB, and JGP designed the concept of the study. PL adjusted the WRF model for LGM applications. KHS performed the regional simulations with the WRF model, analysed the data, and created the figures. XS and GL provided data from the global AWI simulations. PB provided the proxy data. PB and PA contributed to the discussion on the interpretation of proxy data. KHS wrote the first draft of the paper. All authors contributed to discussions and revised the final article.​​​​​​​

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


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


Kim H. Stadelmaier and Joaquim G. Pinto thank the AXA Research Fund for support. Patrick Ludwig is supported by the Helmholtz Climate Initiative REKLIM (regional climate change;, last access: 6 December 2021). Kim H. Stadelmaier, Patrick Ludwig, Xiaoxu Shi, and Gerrit Lohmann thank the German Climate Computing Centre (DKRZ, Hamburg) for providing computing resources. This study is a contribution to the PALEOLINK project (, last access: 6 December 2021) within the PAGES 2k Network; it is also a contribution to the PalMod and PACMEDY projects funded by the BMBF.

Financial support

This research has been supported by the AXA Research Fund (, last access: 6 December 2021​​​​​​​) and the Helmholtz-Gemeinschaft (Climate Initiative REKLIM grant).

The article processing charges for this open-access publication were covered by the Karlsruhe Institute of Technology (KIT).

Review statement

This paper was edited by Alberto Reyes and reviewed by Jef Vandenberghe and one anonymous referee.


Andrieux, E., Bertran, P., Antoine, P., Deschodt, L., Lenoble, A., Coutard, S., Ajas, A., Borderie, Q., Coutard, J.-P., Didierjean, F., Dousteyssier, B., Ferrier, C., Gardère, P., Gé, T., Liard, M., Locht, J.-L., Naton, H.-G., Rué, M., Sitzia, L., Van Vliet-Lanoe, B., and Vernet, G.: Database of pleistocene periglacial featuresin France: description of the online version, Quaternaire, 27, 329–339,, 2016a (data available at:, last access: 6 December 2021)​​​​​​​. a

Andrieux, E., Bertran, P., and Saito, K.: Spatial analysis of the French Pleistocene permafrost by a GIS database, Permafrost. Periglac., 27, 17–30,, 2016b. a, b, c, d, e, f

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,, 2018. a, b, c

Annan, J. D. and Hargreaves, J. C.: A new global reconstruction of temperature changes at the Last Glacial Maximum, Clim. Past, 9, 367–376,, 2013. a, b

Antoine, P., Catt, J., Lautridou, J.-P., and Sommé, J.: The loess and coversands of northern France and southern England, J. Quaternary Sci., 18, 309–318,, 2003. a

Antoine, P., Rousseau, D.-D., Moine, O., Kunesch, S., Hatté, C., Lang, A., Tissoux, H., and Zöller, L.: Rapid and cyclic aeolian deposition during the Last Glacial in European loess: a high-resolution record from Nussloch, Germany, Quaternary Sci. Rev., 28, 2955–2973,, 2009. a

Antoine, P., Goval, E., Jamet, G., Coutard, S., Moine, O., Hérisson, D., Auguste, P., Guérin, G., Lagroix, F., Schmidt, E., Robert, V., Debenham, N., Meszner, S., and Bahain, J.-J.: Les séquences loessiques pléistocène supérieur d'Havrincourt (Pas-de-Calais, France): stratigraphie, paléoenvironnements, géochronologie et occupations paléolithiques, Quaternaire, 25, 321–368,, 2014. a

Barosh, P. J.: Frostquakes in New England, Eng. Geol., 56, 389–394, 2000. a

Bartlein, P. J., Harrison, S. P., Brewer, S., Connor, S., Davis, B. A. S., Gajewski, K., Guiot, J., Harrison-Prentice, T. I., Henderson, A., Peyron, O., Prentice, I. C., Scholze, M., Seppä, H., Shuman, B., Sugita, S., Thompson, R. S., Viau, A. E., Williams, J., and Wu, H.: Pollen-based continental climate reconstructions at 6 and 21 ka: a global synthesis, Clim. Dynam., 37, 775–802,, 2011. a, b, c

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,, 2014. a

Bertran, P., Andrieux, E., Bateman, M., Font, M., Manchuel, K., and Sicilia, D.: Features caused by ground ice growth and decay in Late Pleistocene fluvial deposits, Paris Basin, France, Geomorphology, 310, 84–101,, 2018. a

Bockheim, J. G., Kurz, M. D., Soule, S. A., and Burke, A.: Genesis of active sand-filled polygons in lower and central Beacon Valley, Antarctica, Permafrost. Periglac., 20, 295–308, 2009. a

Braconnot, P., Harrison, S., Kageyama, M., Bartlein, P., Masson-Delmotte, V., Abe-Ouchi, A., Otto-Bliesner, B., and Zhao, Y.: Evaluation of climate models using palaeoclimate data, Nat. Clim. Change, 2, 417–424,, 2012. a, b, c

Brierley, C. M., Zhao, A., Harrison, S. P., Braconnot, P., Williams, C. J. R., Thornalley, D. J. R., Shi, X., Peterschmitt, J.-Y., Ohgaito, R., Kaufman, D. S., Kageyama, M., Hargreaves, J. C., Erb, M. P., Emile-Geay, J., D'Agostino, R., Chandan, D., Carré, M., Bartlein, P. J., Zheng, W., Zhang, Z., Zhang, Q., Yang, H., Volodin, E. M., Tomas, R. A., Routson, C., Peltier, W. R., Otto-Bliesner, B., Morozova, P. A., McKay, N. P., Lohmann, G., Legrande, A. N., Guo, C., Cao, J., Brady, E., Annan, J. D., and Abe-Ouchi, A.: Large-scale features and evaluation of the PMIP4-CMIP6 midHolocene simulations, Clim. Past, 16, 1847–1872,, 2020. a

Burn, C. R.: Implications for palaeoenvironmental reconstruction of recent ice-wedge development at Mayo, Yukon Territory, Permafrost. Periglac., 1, 3–14, 1990. a

Chen, F. and Dudhia, J.: Coupling an Advanced Land Surface–Hydrology Model with the Penn State–NCAR MM5 Modeling System. Part I: Model Implementation and Sensitivity, Mon. Weather Rev., 129, 569–585,<0569:CAALSH>2.0.CO;2, 2001. a

Clark, P. U., Dyke, A. S., Shakun, J. D., Carlson, A. E., Clark, J., Wohlfarth, B., Mitrovica, J. X., Hostetler, S. W., and McCabe, A. M.: The Last Glacial Maximum, Science, 325, 710–714,, 2009. a, b, c

Cleator, S. F., Harrison, S. P., Nichols, N. K., Prentice, I. C., and Roulstone, I.: A new multivariable benchmark for Last Glacial Maximum climate simulations, Clim. Past, 16, 699–712,, 2020. a

CLIMAP Project Members: The last interglacial ocean, Quaternary Res., 2, 123–224,, 1984. a

Dietrich, S. and Seelos, K.: The reconstruction of easterly wind directions for the Eifel region (Central Europe) during the period 40.3–12.9 ka BP, Clim. Past, 6, 145–154,, 2010. a

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958,, 2016. a

FESOM: The AWI Climate Model (AWI-CM), FESOM, available at:, last access: 6 December 2021. a

Flückiger, J., Knutti, R., White, J. W. C., and Renssen, H.: Modeled seasonality of glacial abrupt climate events, Clim. Dynam., 31, 633–645,, 2008. a

Fortier, D. and Allard, M.: Frost-cracking conditions, Bylot Island, eastern Canadian Arctic archipelago, Permafrost. Periglac., 16, 145–161,, 2005. a

Frauenfeld, O. W., Zhang, T., and Mccreight, J. L.: Northern Hemisphere freezing/thawing index variations over the twentieth century, Int. J. Climatol., 27, 47–63,, 2007. a

Harris, C., Arenson, L. U., Christiansen, H. H., Etzelmüller, B., Frauenfelder, R., Gruber, S., Haeberli, W., Hauck, C., Hölzle, M., Humlum, O., Isaksen, K., Kääb, A., Kern-Lütschg, M. A., Lehning, M., Matsuoka, N., Murton, J. B., Nötzli, J., Phillips, M., Ross, N., Seppälä, M., Springman, S. M., and Vonder Mühll, D.: Permafrost and climate in Europe: Monitoring and modelling thermal, geomorphological and geotechnical responses, Earth-Sci. Rev., 92, 117–171,, 2009. a

Harrison, S. P., Bartlein, P. J., Izumi, K., Li, G., Annan, J., Hargreaves, J., P., B., and Kageyama, M.: Evaluation of CMIP5 palaeo-simulations to improve climate projections, Nat. Clim. Change, 5, 735–743,, 2015. a

Harry, D. G. and Gozdzik, J. S.: Ice wedges: Growth, thaw transformation, and palaeoenvironmental significance, J. Quaternary Sci., 3, 39–55,, 1988. a

Hatté, C., Fontugne, M., Rousseau, D.-D., Antoine, P., Zöller, L., Laborde, N. T., and Bentaleb, I.: δ13C variations of loess organic matter as a record of the vegetation response to climatic changes during the Weichselian, Geology, 26, 583–586,<0583:CVOLOM>2.3.CO;2, 1998. a

Hatté, C., Antoine, P., Fontugne, M., Lang, A., Rousseau, D.-D., and Zöller, L.: δ13C of Loess Organic Matter as a Potential Proxy for Paleoprecipitation, Quaternary Res., 55, 33–38,, 2001. a

Hong, S.-Y., Dudhia, J., and Chen, S.-H.: A Revised Approach to Ice Microphysical Processes for the Bulk Parameterization of Clouds and Precipitation, Mon. Weather Rev., 132, 103–120,<0103:ARATIM>2.0.CO;2, 2004. a

Hong, S.-Y., Noh, Y., and Dudhia, J.: A New Vertical Diffusion Package with an Explicit Treatment of Entrainment Processes, Mon. Weather Rev., 134, 2318–2341,, 2006. a

Hughes, A. L. C., Gyllencreutz, R., Lohne, O. S., Mangerud, J., and Svendsen, J. I.: The last Eurasian ice sheets – a chronological database and time-slice reconstruction, DATED-1, Boreas, 45, 1–45,, 2015. a

Huijzer, A. and Isarin, R.: 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,, 1997. a

Huijzer, A. S., Vandenberghe, J., van Huissteden, J., and Isarin, R. F. B​​​​​​​: Paleo-periglacial phenomena in Northwestern Europe, Version 1, Table PERIGLAC, Boulder, Colorado USA, NSIDC: National Snow and Ice Data Center [data set],, 1998. a

Iacono, M. J., Delamere, J. S., Mlawer, E. J., Shephard, M. W., Clough, S. A., and Collins, W. D.: Radiative forcing by long-lived greenhouse gases: Calculations with the AER radiative transfer models, J. Geophys. Res.-Atmos., 113, D13103,, 2008. a

International Research Institute for Climate and Society, Columbia Climate School (IRI): CLIMAP Last Glacial Maximum, IRI [data set], available at: (last access: 6 December 2021), 2015. a

IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, UK and New York, NY, USA, 2013. a, b

Intergovernmental Panel on Climate Change (IPCC): IPCC special report on the ocean and cryosphere in a changing climate, available at: (last access: 6 December 2021), 2019. a, b, c, d

Isarin, R. F. B., Renssen, H., and Vandenberghe, J.: The impact of the North Atlantic Ocean on the Younger Dryas climate in northwestern and central Europe, J. Quaternary Sci., 13, 447–453,<447::AID-JQS402>3.0.CO;2-B, 1998. a, b

Jiménez, P. A., Dudhia, J., González-Rouco, J. F., Navarro, J., Montávez, J. P., and García-Bustamante, E.: A Revised Scheme for the WRF Surface Layer Formulation, Mon. Weather Rev., 140, 898–918,, 2012. a

Jungclaus, J., Giorgetta, M., Reick, C., Legutke, S., Brovkin, V., Crueger, T., Esch, M., Fieg, K., Fischer, N., Glushak, K., Gayler, V., Haak, H., Hollweg, H.-D., Kinne, S., Kornblueh, L., Matei, D., Mauritsen, T., Mikolajewicz, U., Müller, W., Notz, D., Pohlmann, T., Raddatz, T., Rast, S., Roeckner, E., Salzmann, M., Schmidt, H., Schnur, R., Segschneider, J., Six, K., Stockhause, M., Wegner, J., Widmann, H., Wieners, K.-H., Claussen, M., Marotzke, J., and Stevens, B.: CMIP5 simulations of the Max Planck Institute for Meteorology (MPI-M) based on the MPI-ESM-P model: The lgm experiment, served by ESGF, WCRP [data set],, 2012a. a

Jungclaus, J., Giorgetta, M., Reick, C., Legutke, S., Brovkin, V., Crueger, T., Esch, M., Fieg, K., Fischer, N., Glushak, K., Gayler, V., Haak, H., Hollweg, H.-D., Kinne, S., Kornblueh, L., Matei, D., Mauritsen, T., Mikolajewicz, U., Müller, W., Notz, D., Pohlmann, T., Raddatz, T., Rast, S., Roeckner, E., Salzmann, M., Schmidt, H., Schnur, R., Segschneider, J., Six, K., Stockhause, M., Wegner, J., Widmann, H., Wieners, K.-H., Claussen, M., Marotzke, J., and Stevens, B.: CMIP5 simulations of the Max Planck Institute for Meteorology (MPI-M) based on the MPI-ESM-P model: The piControl experiment, served by ESGF, WCRP [data set],, 2012b. a

Jungclaus, J. H., Fischer, N., Haak, H., Lohmann, K., Marotzke, J., Matei, D., Mikolajewicz, U., Notz, D., and von Storch, J. S.: Characteristics of the ocean simulations in the Max Planck Institute Ocean Model (MPIOM) the ocean component of the MPI-Earth system model, J. Adv. Model. Earth Sy., 5, 422–446,, 2013. a

Justino, F. and Peltier, W. R.: The glacial North Atlantic Oscillation, Geophys. Res. Lett., 32, L21803,, 2005. a

Kageyama, M., Albani, S., Braconnot, P., Harrison, S. P., Hopcroft, P. O., Ivanovic, R. F., Lambert, F., Marti, O., Peltier, W. R., Peterschmitt, J.-Y., Roche, D. M., Tarasov, L., Zhang, X., Brady, E. C., Haywood, A. M., LeGrande, A. N., Lunt, D. J., Mahowald, N. M., Mikolajewicz, U., Nisancioglu, K. H., Otto-Bliesner, B. L., Renssen, H., Tomas, R. A., Zhang, Q., Abe-Ouchi, A., Bartlein, P. J., Cao, J., Li, Q., Lohmann, G., Ohgaito, R., Shi, X., Volodin, E., Yoshida, K., Zhang, X., and Zheng, W.: The PMIP4 contribution to CMIP6 – Part 4: Scientific objectives and experimental design of the PMIP4-CMIP6 Last Glacial Maximum experiments and PMIP4 sensitivity experiments, Geosci. Model Dev., 10, 4035–4055,, 2017. a

Kageyama, M., Harrison, S. P., Kapsch, M.-L., Lofverstrom, M., Lora, J. M., Mikolajewicz, U., Sherriff-Tadano, S., Vadsaria, T., Abe-Ouchi, A., Bouttes, N., Chandan, D., Gregoire, L. J., Ivanovic, R. F., Izumi, K., LeGrande, A. N., Lhardy, F., Lohmann, G., Morozova, P. A., Ohgaito, R., Paul, A., Peltier, W. R., Poulsen, C. J., Quiquet, A., Roche, D. M., Shi, X., Tierney, J. E., Valdes, P. J., Volodin, E., and Zhu, J.: The PMIP4 Last Glacial Maximum experiments: preliminary results and comparison with the PMIP3 simulations, Clim. Past, 17, 1065–1089,, 2021. a, b, c

Kain, J. S.: The Kain–Fritsch Convective Parameterization: An Update, J. Appl. Meteorol., 43, 170–181,<0170:TKCPAU>2.0.CO;2, 2004. a

Keeble, J., Hassler, B., Banerjee, A., Checa-Garcia, R., Chiodo, G., Davis, S., Eyring, V., Griffiths, P. T., Morgenstern, O., Nowack, P., Zeng, G., Zhang, J., Bodeker, G., Burrows, S., Cameron-Smith, P., Cugnet, D., Danek, C., Deushi, M., Horowitz, L. W., Kubin, A., Li, L., Lohmann, G., Michou, M., Mills, M. J., Nabat, P., Olivié, D., Park, S., Seland, Ø., Stoll, J., Wieners, K.-H., and Wu, T.: Evaluating stratospheric ozone and water vapour changes in CMIP6 models from 1850 to 2100, Atmos. Chem. Phys., 21, 5015–5061,, 2021. a

Kitover, D. C., van Balen, R. T., Roche, D. M., Vandenberghe, J., and Renssen, H.: New Estimates of Permafrost Evolution during the Last 21 k Years in Eurasia using Numerical Modelling, Permafrost. Periglac., 24, 286–303,, 2013. a

Kokelj, S., Lantz, T., Wolfe, S. A., Kanigan, J., Morse, P., Coutts, R., Molina-Giraldo, N., and Burn, C. R.: Distribution and activity of ice wedges across the forest-tundra transition, western Arctic Canada, J. Geophys. Res.-Earth, 119, 2032–2047, 2014. a, b

Krauß, L., Zens, J., Zeeden, C., Schulte, P., Eckmeier, E., and Lehmkuhl, F.: A Multi-Proxy Analysis of two Loess-Paleosol Sequences in the Northern Harz Foreland, Germany, Palaeogeography, Palaeoclimatology, Palaeoecology, 461, 401–417,, 2016. a, b

Lambeck, K., Yokoyama, Y., and Purcell, T.: Into and out of the Last Glacial Maximum: sea-level change during Oxygen Isotope Stages 3 and 2, Quaternary Sci. Rev., 21, 343–360,, ePILOG, 2002. a

Lambeck, K., Rouby, H., Purcell, A., Sun, Y., and Sambridge, M.: Sea level and global ice volumes from the Last Glacial Maximum to the Holocene, P. Natl. Acad. Sci. USA, 111, 15296–15303, 2014. a

Lehmkuhl, F., Zens, J., Krauß, L., Schulte, P., and Kels, H.: Loess-paleosol sequences at the northern European loess belt in Germany: Distribution, geomorphology and stratigraphy, Quaternary Sci. Rev., 153, 11–30,, 2016. a

Levavasseur, G., Vrac, M., Roche, D. M., Paillard, D., Martin, A., and Vandenberghe, J.: Present and LGM permafrost from climate simulations: contribution of statistical downscaling, Clim. Past, 7, 1225–1246,, 2011. a, b

Levy, J. S., Head, J. W., and Marchant, D. R.: The role of thermal contraction crack polygons in cold-desert fluvial systems, Antarct. Sci., 20, 565–579,, 2008. a

Li, C. and Battisti, D. S.: Reduced Atlantic Storminess during Last Glacial Maximum: Evidence from a Coupled Climate Model, J. Climate, 21, 3561–3579,, 2008. a, b

Liu, Y. and Jiang, D.: Last glacial maximum permafrost in China from CMIP5 simulations, Palaeogeography, Palaeoclimatology, Palaeoecology, 447, 12–21,, 2016a. a, b

Liu, Y. and Jiang, D.: Mid-Holocene permafrost: Results from CMIP5 simulations, J. Geophys. Res.-Atmos., 121, 221–240,, 2016b. a, b, c

Lohmann, G., Butzin, M., Eissner, N., Shi, X., and Stepanek, C.: Abrupt Climate and Weather Changes Across Time Scales, Paleoceanography and Paleoclimatology, 35, e2019PA003782,, 2020. a, b

Ludwig, P., Schaffernicht, E. J., Shao, Y., and Pinto, J. G.: Regional atmospheric circulation over Europe during the Last Glacial Maximum and its links to precipitation, J. Geophys. Res.-Atmos., 121, 2130–2145,, 2016. a, b, c, d, e

Ludwig, P., Pinto, J. G., Raible, C. C., and Shao, Y.: Impacts of surface boundary conditions on regional climate model simulations of European climate during the Last Glacial Maximum, Geophys. Res. Lett., 44, 5086–5095,, 2017. a, b, c, d, e

Ludwig, P., Gómez-Navarro, J. J., Pinto, J. G., Raible, C. C., Wagner, S., and Zorita, E.: Perspectives of regional paleoclimate modeling, Ann. NY Acad. Sci., 1436, 54–69,, 2019. a, b

Löfverström, M., Caballero, R., Nilsson, J., and Kleman, J.: Evolution of the large-scale atmospheric circulation in response to changing ice sheets over the last glacial cycle, Clim. Past, 10, 1453–1471,, 2014. a, b, c

Mackay, J. R. and Burn, C. R.: The first 20 years (1978–1979 to 1998–1999) of active-layer development, Illisarvik experimental drained lake site, western Arctic coast, Canada1, Can. J. Earth Sci., 39, 1657–1674, 2002. a

MARGO Project Members: Constraints on the magnitude and patterns of ocean cooling at the Last Glacial Maximum, Nat. Geosci., 2, 127–132,, 2009. a

Marsland, S. J., Haak, H., Jungclaus, J. H., Latif, M., and Röske, F.: The Max-Planck-Institute global ocean/sea ice model with orthogonal curvilinear coordinates, Ocean Model., 5, 91–127,, 2003. a

Matsuoka, N., Christiansen, H. H., and Watanabe, T.: Ice-wedge polygon dynamics in Svalbard: Lessons from a decade of automated multi-sensor monitoring, Permafrost. Periglac., 29, 210–227,, 2018. a, b, c, d, e, f, g, h

Merz, N., Raible, C. C., and Woollings, T.: North Atlantic Eddy-Driven Jet in Interglacial and Glacial Winter Climates, J. Climate, 28, 3977–3997,, 2015. a, b

Meszner, S., Kreutzer, S., Fuchs, M., and Faust, D.: Late Pleistocene landscape dynamics in Saxony, Germany: Paleoenvironmental reconstruction using loess-paleosol sequences, Quatern. Int., 296, 94–107,, 2013. a

Mix, A. C., Bard, E., and Schneider, R.: Environmental processes of the ice age: land, oceans, glaciers (EPILOG), Quaternary Sci. Rev., 20, 627–657,, 2001. a

Monnin, E., Indermühle, A., Dällenbach, A., Flückiger, J., Stauffer, B., Stocker, T. F., Raynaud, D., and Barnola, J.-M.: Atmospheric CO2 Concentrations over the Last Glacial Termination, Science, 291, 112–114,, 2001. a

Murton, J. B.: Permafrost and Periglacial Features. Ice wedges and ice wedge casts, in: Encyclopedia of Quaternary Science, edited by: Elias, S. A., vol. 2, Elsevier, Amsterdam, 436–451, 2013. a

Murton, J. B., Worsley, P., and Gozdzik, J.: Sand veins and wedges in cold aeolian environments, Quaternary Sci. Rev., 19, 899–922,, 2000. a

Nelson, F. E. and Outcalt, S. I.: A Computational Method for Prediction and Regionalization of Permafrost, Arctic Alpine Res., 19, 279–288, DOI: 10.1080/00040851.1987.12002602, 1987. a, b

Niu, G.-Y., Yang, Z.-L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari, M., and Xia, Y.: The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements, J. Geophys. Res.-Atmos., 116, D12109,, 2011. a

Okkonen, J., Neupauer, R., Kozlovskaya, E., Afonin, N., Moisio, K., Taewook, K., and Muurinen, E.: Frost Quakes: Crack Formation by Thermal Stress, J. Geophys. Res.-Earth, 125, e2020JF005616,, 2020. a

Pausata, F. S. R., Li, C., Wettstein, J. J., Kageyama, M., and Nisancioglu, K. H.: The key role of topography in altering North Atlantic atmospheric circulation during the last glacial period, Clim. Past, 7, 1089–1101,, 2011. a

Peltier, W. R., Argus, D. F., and Drummond, R.: Space geodesy constrains ice age terminal deglaciation: The global ICE-6G_C (VM5a) model, J. Geophys. Res.-Sol. Ea., 120, 450–487,, 2015. a, b

Péwé, T. L.: Sand-wedge polygons (tesselations) in the McMurdo Sound region, Antarctica; a progress report, Am. J. Sci., 257, 545–552, 1959. a

Péwé, T. L.: Paleoclimatic significance of fossil ice wedges, Biuletyn Peryglacjalny, 15, 65–73, 1966. a, b

Pfahl, S., O'Gorman, P. A., and Singh, M. S.: Extratropical Cyclones in Idealized Simulations of Changed Climates, J. Climate, 28, 9373–9392,, 2015. a

Pinto, J. G. and Ludwig, P.: Extratropical cyclones over the North Atlantic and western Europe during the Last Glacial Maximum and implications for proxy interpretation, Clim. Past, 16, 611–626,, 2020. a, b

Prentice, I. C. and Harrison, S. P.: Ecosystem effects of CO2 concentration: evidence from past climates, Clim. Past, 5, 297–307,, 2009. a

Prospero, J. M., Ginoux, P., Torres, O., Nicholson, S. E., and Gill, T. E.: Environmental characterization of global sources of atmoshperic soil dust identified with the NIMBUS 7 total ozone mapping spectrometer (TOMS) absorbing aerosol product, Rev. Geophys., 40, 1002,, 2002. a

Raible, C. C., Pinto, J. G., Ludwig, P., and Messmer, M.: A review of past changes in extratropical cyclones in the northern hemisphere and what can be learned for the future, WIREs Clim. Change, 12, E680,, 2021. a, b

Ray, N. and Adams, J. M.: A GIS-based vegetation map of the world at the Last Glacial Maximum (25,000–15,000 BP), Internet Archaeology, 11, 2001. a

Renssen, H., Kasse, C., Vandenberghe, J., and Lorenz, S. J.: Weichselian Late Pleniglacial surface winds over northwest and central Europe: a model–data comparison, J. Quaternary Sci., 22, 281–293,, 2007. a

Romanovskij, N.: Regularities in formation of frost fissures and development of frost fissure polygons, Biuletyn Peryglacjalny, 23, 237–277, 1973. a

Royer, A., Picard, G., Vargel, C., Langlois, A., Gouttevin, I., and Dumont, M.: Improved Simulation of Arctic Circumpolar Land Area Snow Properties and Soil Temperatures, Front. Earth Sci., 9, 515,, 2021. a

Römer, W., Lehmkuhl, F., and Sirocko, F.: Late Pleistocene aeolian dust provenances and wind direction changes reconstructed by heavy mineral analysis of the sediments of the Dehner dry maar (Eifel, Germany), Global Planet. Change, 147, 25–39,, 2016. a

Saito, K., Sueyoshi, T., Marchenko, S., Romanovsky, V., Otto-Bliesner, B., Walsh, J., Bigelow, N., Hendricks, A., and Yoshikawa, K.: LGM permafrost distribution: how well can the latest PMIP multi-model ensembles perform reconstruction?, Clim. Past, 9, 1697–1714,, 2013. a, b

Schaffernicht, E. J., Ludwig, P., and Shao, Y.: Linkage between dust cycle and loess of the Last Glacial Maximum in Europe, Atmos. Chem. Phys., 20, 4969–4986,, 2020. a, b

Schuur, E. A. G., McGuire, A. D., Schädel, C., Grosse, G., Harden, J. W., Hayes, D. J., Hugelius, G., D., K. C., Kuhry, P., Lawrence, D. M., Natali, S. M., Olefeldt, D., Romanovsky, V. E., Schaefer, K., Turetsky, M. R., Treat, C. C., and Vonk, J. E.: Climate change and the permafrost carbon feedback, Nature, 520, 171–179,, 2015. a, b

Schwan, J.: The origin of horizontal alternating bedding in weichselian aeolian sands in Northwestern Europe, Sediment. Geol., 49, 73–108,, 1986. a

Schwan, J.: The structure and genesis of Weichselian to early hologene aeolian sand sheets in western Europe, Sediment. Geol., 55, 197–232,, 1988. a

Seguinot, J., Ivy-Ochs, S., Jouvet, G., Huss, M., Funk, M., and Preusser, F.: Modelling last glacial cycle ice dynamics in the Alps, The Cryosphere, 12, 3265–3285,, 2018. a

Sidorenko, D., Rackow, T., Semmler, T., Barbi, D., Danilov, S., Dethloff, K., Dorn, W., Fieg, K., Goessling, H. F., Handorf, D., Harig, S., Hiller, W., Juricke, S., Losch, M., Schröter, J., Sein, D. V., and Wang, Q.: Towards multi-resolution global climate modeling with ECHAM6–FESO M. Part I: model formulation and mean climate, Clim. Dynam., 44, 757–780,, 2015. a

Sima, A., Rousseau, D.-D., Kageyama, M., Ramstein, G., Schulz, M., Balkanski, Y., Antoine, P., Dulac, F., and Hatté, C.: Imprint of North-Atlantic abrupt climate changes on western European loess deposits as viewed in a dust emission model, Quaternary Sci. Rev., 28, 2851–2866,, 2009. a

Skamarock, W., Klemp, J., Dudhia, J., Gill, D., Barker, D., Duda, M., Huang, X.-Y., Wang, W., and Powers, J.: A description of the advanced research WRF version 3, NCAR Tech. Note NCAR/TN-475+STR, 113,, 2008. a

Slater, A. G. and Lawrence, D. M.: Diagnosing Present and Future Permafrost from Climate Models, J. Climate, 26, 5608–5623,, 2013. a, b

Smerdon, J. E.: Climate models as a test bed for climate reconstruction methods: pseudoproxy experiments, WIREs Clim. Change, 3, 63–77,, 2012. a

Smith, M. and Riseborough, D.: Climate and the limits of permafrost: a zonal analysis, Permafrost. Periglac., 13, 1–15, 2002. a, b

Stendel, M. and Christensen, J. H.: Impact of global warming on permafrost conditions in a coupled GCM, Geophys. Res. Lett., 29, 1632,, 2002. a, b

Stevens, B., Giorgetta, M., Esch, M., Mauritsen, T., Crueger, T., Rast, S., Salzmann, M., Schmidt, H., Bader, J., Block, K., Brokopf, R., Fast, I., Kinne, S., Kornblueh, L., Lohmann, U., Pincus, R., Reichler, T., and Roeckner, E.: Atmospheric component of the MPI-M Earth System Model: ECHAM6, J. Adv. Model. Earth Sy., 5, 146–172,, 2013. a

Stevens, T., Sechi, D., Bradák, B., Orbe, R., Baykal, Y., Cossu, G., Tziavaras, C., Andreucci, S., and Pascucci, V.: Abrupt last glacial dust fall over southeast England associated with dynamics of the British-Irish ice sheet, Quaternary Sci. Rev., 250, 106641,, 2020. a, b

Tarasov, L. and Peltier, R., W.: Greenland glacial history, borehole constraints, and Eemian extent, J. Geophys. Res.-Sol. Ea., 108, 2143,, 2003. a

Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An Overview of CMIP5 and the Experiment Design, B. Am. Meteorol. Soc., 93, 485–498,, 2012. a

Tewari, M., Chen, F., Wang, W., Dudhia, J., LeMone, M. A., Mitchell, K., Ek, M., Gayno, G., Wegiel, J., and Cuenca, R. H.: Implementation and verification of the unified NOAH land surface model in the WRF model, in: 20th conference on weather analysis and forecasting/16th conference on numerical weather prediction, Seattle, WA, 10 to 16 January 2004, available at: (last access: 6 December 2021)​​​​​​​, 11–15, 2004. a, b

Throop, J., Lewkowicz, A. G., and Smith, S. L.: Climate and ground temperature relations at sites across the continuous and discontinuous permafrost zones, northern Canada, Can. J. Earth Sci., 49, 865–876, 2012. a

Ullman, D. J., LeGrande, A. N., Carlson, A. E., Anslow, F. S., and Licciardi, J. M.: Assessing the impact of Laurentide Ice Sheet topography on glacial climate, Clim. Past, 10, 487–507,, 2014. a

van Everdingen, R. O.: Multi-language glossary of permafrost and related ground-ice terms, National Snow and Ice Data Center/World Data Center for Glaciology, Boulder, CO, 2005. a

Vandenberghe, J.: Some periglacial phenomena and their stratigraphical position in Weichselian deposits in the Netherlands, Polarforschung, 53, 97–107, 1983. a

Vandenberghe, J., Renssen, H., Roche, D. M., Goosse, H., Velichko, A. A., Gorbunov, A., and Levavasseur, G.: Eurasian permafrost instability constrained by reduced sea-ice cover, Quaternary Sci. Rev., 34, 16–23,, 2012.  a, b

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,, 2014. a, b, c, d

Vandenberghe, J., French, H. M., Jin, H., Wang, X., Yi, S., and He, R.: The extent of permafrost during the Last Permafrost Maximum (LPM) on the Ordos Plateau, north China, Quaternary Sci. Rev., 214, 87–97,, 2019. a

Wang, Q., Danilov, S., Sidorenko, D., Timmermann, R., Wekerle, C., Wang, X., Jung, T., and Schröter, J.: The Finite Element Sea Ice-Ocean Model (FESOM) v.1.4: formulation of an ocean general circulation model, Geosci. Model Dev., 7, 663–693,, 2014. a

Wang, T., Liu, Y., and Huang, W.: Last Glacial Maximum Sea Surface Temperatures: A Model-Data Comparison, Atmospheric and Oceanic Science Letters, 6, 233–239,, 2013. a

Washburn, A. L.: Frost cracking in a middle latitude climate, Biuletyn Peryglacjalny, 12, 175–189, 1963. a

Washburn, A. L.: Geocryology: a survey of periglacial processes and environments, Wiley, Hoboken, NJ, USA, 1979. a

Woillez, M.-N., Kageyama, M., Krinner, G., de Noblet-Ducoudré, N., Viovy, N., and Mancip, M.: Impact of CO2 and climate on the Last Glacial Maximum vegetation: results from the ORCHIDEE/IPSL models, Clim. Past, 7, 557–577,, 2011. a

Wolfe, S. A., Morse, P. D., Neudorf, C. M., Kokelj, S. V., Lian, O. B., and O'Neill, H. B.: Contemporary sand wedge development in seasonally frozen ground and paleoenvironmental implications, Geomorphology, 308, 215–229,, 2018. a

Short summary
We use regional climate simulations for the Last Glacial Maximum to reconstruct permafrost and to identify areas of thermal contraction cracking of the ground in western Europe. We find ground cracking, a precondition for the development of permafrost proxies, south of the probable permafrost border, implying that permafrost was not the limiting factor for proxy development. A good agreement with permafrost and climate proxy data is achieved when easterly winds are modelled more frequently.