The effects of late Cenozoic climate change on the global distribution of frost cracking

. Frost cracking is a dominant mechanical weathering phenomenon facilitating the breakdown of bedrock in periglacial regions. Despite recent advances in understanding frost cracking processes, few studies have addressed how global climate change over the late Cenozoic may have impacted spatial variations in frost cracking intensity. In this study, we estimate global changes in frost cracking intensity (FCI) by segregation ice growth. Existing process-based models of FCI are applied in combination with soil thickness data from the Harmonized World Soil Database. Temporal and spatial variations in FCI are predicted using surface temperature changes obtained from ECHAM5 general circulation model simulations conducted for four different paleoclimate time slices. Time slices considered include pre-industrial ( ∼ 1850 CE, PI), mid-Holocene ( ∼ 6 ka, MH), Last Glacial Maximum ( ∼ 21 ka, LGM), and Pliocene ( ∼ 3 Ma, PLIO) times. Results indicate for all paleoclimate time slices that frost cracking was most prevalent (relative to PI times) in the middle- to high-latitude regions, as well as high-elevation lower-latitude areas such the Himalayas, Tibet, the European Alps, the Japanese Alps, the US Rocky Mountains, and the Andes Mountains. The smallest deviations in frost cracking (relative to PI conditions) were observed in the MH simulation, which yielded slightly higher FCI values in most of the areas. In contrast, larger deviations were observed in the simulations of the colder climate (LGM) and warmer climate (PLIO). Our results indicate that the impact of climate change on frost cracking was most severe during the PI–LGM period due to higher differences in temperatures and glaciation at higher latitudes. The PLIO results indicate low FCI in the Andes and higher values of FCI in Greenland and Canada due to the diminished extent of glaciation in the warmer PLIO climate.


94
Soil thickness data was obtained from the re-gridded Harmonized World Soil Database (HWSD) v1.2  which 95 has a 0.05-degree spatial resolution and depths ranging from 0 m to 1 m (Fig. 1)   where, T represents daily subsurface temperature at depth z (m) and time t (days in a year), MAT and Ta represent mean annual 122 surface temperature and half amplitude of annual temperature variation respectively, Py is the period of the sinusoidal cycle (1 123

132
The calculation of subsurface temperatures was discretized into 200 depth intervals from the surface to the maximum depth of 133 20 m. Smaller depth intervals (~1 cm) were used near the surface and large intervals (~20 cm) at greater depths, because the 134 FCI is expected to change most dramatically near the surface and dampen with depth due to thermal diffusion (Andersen et 135 al., 2015). 136

Estimation of Frost Cracking Intensity 137
We applied three different approaches (models) with different levels of complexity to estimate global variations in frost 138 cracking during different past climates (Fig 4; Andersen et al., 2015; Anderson, 1998;Hales and Roering, 2007). The models 139 use predicted ground surface temperatures from each grid cell in the GCM to calculate subsurface temperatures and FCI. We 140 then calculate differences between the FCI from the PI reference simulation and the FCI predicted for the PLIO, LGM and 141 in the models used in our study, which are discussed in detail in sections 3.2.1 -3.2.3. Models 1-3 successively increase in 143 complexity and consider more factors. The approach of Andersen et al., (2015), referred to here as Model 3, is the most recent 144 and complete in its consideration of the processes (e.g. effect of soil-cover on FCI) that are relevant for frost cracking. Given 145 this, we focus our presentation of results in the main text here on Model 3, but for completeness describe below differences of 146 Model 3 from earlier Models (1-2). For brevity, results from the earlier models are presented in the supplementary material. A 147 flowchart illustrating our methods is presented in Fig. 5. Similar to previous studies, the hydrogeological properties of the 148 bedrock (i.e. infiltration, water saturation, porosity and permeability) are ignored in this study. This approach provides a 149 simplified means for estimating the FCI for underlying bedrock at a global scale.    Model 1 represents the simplest approach and applies the method of Anderson (1998). In our application of this model we use 160 a more representative thermal diffusivity value for rocks of 1.5 × 10 −7 2 −1 , because the previous study was specific to 161 granitic bedrock and applied a diffusivity specific to that. Furthermore, the boundary conditions of a low rock surface albedo 162 (≤ 0.1) and presence of a high atmospheric transmissivity (≥ 0.9) on the surface were relaxed, as surface temperatures were 163 used in our study instead of near-surface air temperatures. 164 This implementation supports frost cracking in the bedrock with temperatures between -8 °C and -3 °C (Hallet et al., 1991). 185 In the case if permafrost areas, MAT is always negative, but as sinusoidal T (z, t) is calculated based on MAT and Ta, a positive 186 T (> 0 °C) may occur during warmer days of the year. In addition, Ta is higher for higher latitudes (Fig. 2), which are more 187 prone to frost cracking. 188 The model is described as follows: 189 where FCI (z, t) is the frost cracking intensity at depth z and time t. It is an index for the absolute value of the thermal gradient 192 at that particular depth and time that fulfils the conditions defined above. 193 In equation 5, ´ represents the integrated FCI for a geographic location. More specifically, the FCI is integrated over one

Model 3: Frost cracking intensity as a function of thermal gradients and sediment thickness 198
In the final (most complex) approach used in this study, the effect of an overlying soil layer (Fig. 1)   In Model 3, frost cracking intensity is estimated as a product of the thermal gradient and volume of water available (VW) for 206 segregation ice growth at each depth element, such that: 207    (30 °C to 40 °C) tend to be higher where sediment cover is thicker (e.g. North East Eurasia). 272

Discussion 273
In this section, we synthesize and interpret the results (from Model 3) at both a global and regional scales. For brevity, we limit 274 our discussion regional variations to include Tibet, Europe and South America. For other regional areas of interest to readers, 275 the data used in the following figures is available for download (see acknowledgements). Our presentation of selected regional 276 https://doi.org/10.5194/esurf-2021-78 Preprint. Discussion started: 27 January 2022 c Author(s) 2022. CC BY 4.0 License. areas is followed by the comparison of modeled FCI with published field observations and permafrost extent in the LGM and 277 present day. We also compare the model outcomes of all the three models used in the study. relative to the PI climate simulation on a global scale (Fig. 7), in Tibet (Fig. 8), Europe (Fig. 9) and South America (Fig. 10). 297 Furthermore, the spatial distribution of FCI in various climates has been compared with the glacier mask (Supplement Fig. 3) 298 where continental ice was located for all time-slices (PI, MH, LGM and PLIO). This was done to understand the reasons behind 299 the trend of FCI over time. 300 The differences in FCI between the PI and MH climate simulations are the range of -0.04 °C m to 0.02 °C m. The MH 301 simulation yields s higher FCI values for most regions except for parts of northern Asia, mid-western Europe, mid North 302 America, Alaska, the Andes Mountains and Tibet. These differences may be attributed to the slight changes in MATs in these 303

regions. 304
The differences between PI and LGM FCI values are highest in the high latitudes (Fig. 7) in North America (ΔFCI ≈ 0.16 °C 305 m) and northern Europe (ΔFCI ≈ 0.04 °C m). This is likely due to continental glaciation in these areas (Supplement Fig. 3) 306 leading to low or no frost cracking during LGM. In the mid-high latitudes (~ 50 °N to 70 °N) of Northern Asia, LGM FCI 307 values are higher than in PI FCI values (ΔFCI ≈ -0.06 °C m). This can be attributed to an absence of glacial cover and higher 308 Ta values (Fig. 3)   In summary, the overall comparison of differences between paleo-FCI and PI-FCI indicates a low impact of changing surface 324 temperatures between the PI and MH simulations on frost cracking. This is not surprising given the relatively small 325 climatological differences between the simulations. The differences in FCI between PLIO and PI are more varied, but generally 326 greater. More specifically, warmer regional climates in the Pliocene seem to facilitate frost cracking in higher latitudes, 327 especially in Greenland, when glaciation is absent. The LGM simulation produced the greatest differences in FCI with respect 328 to the PI simulation. These can be attributed to increased glaciation and a much colder climate in higher latitudes, including 329 North America and Europe. High LGM-FCI values were exhibited east of the Andes Mountains in the southern part of South 330 America, possibly due to the absence of glacial cover and high Ta values (~ 20 °C -25 °C) (Fig. 3)

365
The magnitude of PI-MH FCI differences in southwestern South America (Fig. 10)

406
High LGM-FCI values (> ~0.06) correlate with permafrost extent in higher latitudes in Eurasia, Alaska and Tibet and 407 Himalayan regions (Fig. 11). Similarly, a good correlation for FCI > ~-0.06 in the PI simulation and permafrost in the PD 408 simulation is observed for higher latitudes in Asia and North America (Fig 12). However, permafrost extent does not cover 409 Montane and Alpine permafrost in North America near the southern extent of the ice sheet (35 °N to 45 °N latitudes) in the 410 LGM simulation, as suggested by French and Millar (2014), where FCI exhibits a sudden increase ( approx. 0.08 °C m to 0.1 411 °C m). Overall, FCI and permafrost extent correlate reasonably well with our predictions.

5.5.
Inter-comparison of Models 1-3 416 A comparison of the FCI predicted by the three models for the different time slices highlights some key differences (Fig. 6,  417 and supplement Figs 1, 2). The spatial extent of frost cracking in specific time-slices is different in all the three models, which 418 can be accredited to different inputs considered in each model, namely the availability of water for frost cracking by segregation 419 ice growth. For example, the glaciated region (Supplement Fig. 3)