Heterogeneous changes of soil microclimate in high mountains and glacier forelands

doi:10.21203/rs.3.rs-2017904/v1

Download PDF

Article

Heterogeneous changes of soil microclimate in high mountains and glacier forelands

https://doi.org/10.21203/rs.3.rs-2017904/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 31 Aug, 2023

Read the published version in Nature Communications →

Version 1

posted

You are reading this latest preprint version

Landscapes nearby glaciers are disproportionally affected by climate change, still we lack the information on microclimate variation that is required to understand impacts of climate change on these ecosystems and their biodiversity. Here we use near-subsurface soil temperatures in 175 stations from polar, equatorial and alpine glacier forelands to reconstruct temperatures at high resolution, assess spatial differences in microclimate change from 2001 to 2020, and estimate whether microclimate heterogeneity might buffer the severity of warming impacts on organisms. Temporal changes in microclimate are tightly linked to broad-scale trends, but the rate of global warming showed spatial heterogeneity, with faster warming nearby glaciers and during the warm season, and an extension of the snow-free season. Still, the fine-scale spatial variability of microclimate is one-to-ten times larger than the temporal change experienced, indicating the potential for microclimate to buffer climate change, possibly allowing organism to withstand, at least temporarily, the effects of warming.

soil temperature

glacier foreland

microclimate

microclimatic buffering

global changes

Mountain ecosystems provide multiple goods and services to humankind and act as fundamental regulators of regional climate and hydrology^1,2,3. The topographic and climatic heterogeneity of mountain areas, as well as their geological history, deeply influence several biological processes (i.e. adaptation, speciation, dispersal, persistence, and extinction⁴); as a result, mountain ecosystems are biodiversity hotspots with unique levels of endemism, adaptations and lifeforms⁵. Despite mountains cover only one-tenth of the Earth continental surface (excluding Antarctica)⁶, one-quarter of all terrestrial species live there⁷. However, ongoing climatic changes are causing unprecedented modifications of mountain systems^2,3. At the highest elevations, glaciers are losing mass, and the pace of glacier retreat has been globally accelerating during the past decades⁸. This dramatic glacier shrinkage has multiple impacts on all biotic and abiotic components of ecosystems^3,9,10,11,12. Recently deglaciated lands undergo rapid geomorphological transformations with loose sediments from the early-successional stages rapidly developing into structured soils^9,13,14; in turn, this facilitates the colonization of the foreland by multiple lifeforms^12,15. However, climatic variations affect the rates of change of these ecosystems. For instance, temperature influences the rate of rock weathering^16,17 and warmer areas experience faster soil development¹⁴. High temperatures can also affect temporal dynamics of communities, by increasing colonization success by termophilic species and favouring evolution towards more complex community structures^18,19,20, and influence carbon fluxes between soil, vegetation and the atmosphere^21,22.

In areas with complex topography, regional climate (i.e. macroclimate) interacts with topography, potentially resulting in local temperatures decoupled from the regional average²³. Microclimate can be defined as the fine-scale spatial and temporal offsets of the local climate from the macroclimate²⁴. The decoupling between micro- and macroclimate is particularly pronounced near, and below, soil surface^25,26, where microclimate best represents the set of climatic conditions actually experienced by organisms. In mountain areas, topographic elements (i.e. elevation, aspect, slope and topographic shading) locally regulate incoming solar radiation, evapotranspiration, wind speed, cold air drainage, and snow accumulation and melt at fine spatial scales, generating complex patterns of local climatic conditions^23,27,28. Along glacier forelands, additional factors influence local climate, such as vegetation cover and height, distance from the ice mass, and soil texture, creating heterogeneous microhabitats inhabited by different biotic assemblages²⁹. Snow cover is a further key driver of the functioning of mountain ecosystems³⁰, affecting biogeochemical and hydrological processes, and controlling the life cycle of many organisms by determining the duration of their growing / activity season^30,31,32, with potential impacts on ecosystem productivity³¹.

At fine spatial scales, spatial variability of both local temperature and snow can be strong, creating a mosaic of nearby micro-habitats that host different communities^31,33. Microclimate differences between nearby areas might at least temporarily buffer the severity of warming impacts on populations. Microclimate buffering is the dampening of macro-climatic fluctuations due to local conditions (e.g. topography and vegetation cover), such that these fluctuations still exist at the microclimatic scale, but have lower intensity and a more limited effect on organisms³⁴. Microclimate heterogeneity can buffer the impacts of macroclimate change if individuals are able to move between neighbouring sites characterized by microclimatic differences^35,36,37,38. Detailed information on microclimate and snow cover is thus pivotal to understand the consequences of climate change on organisms living in mountain ecosystems and the potential buffering effect favoured by microclimatic heterogeneity, still this requires global scale, high-resolution analyses that were so far lacking.

Here, we use a unique dataset of near-subsurface soil temperatures collected in 175 stations from polar, equatorial and alpine glacier forelands to produce a high-resolution, global reconstruction of monthly average soil temperatures during the snow-free season in high mountains and proglacial environments (Fig. 1). To combine the accuracy of empirically-calibrated relationships with the transferability of mechanistic models, we implemented a correlative hybrid approach based on the mechanistic understanding of the main drivers of microclimate^39,40. Terms were introduced in the modelling framework to account for both the horizontal (elevation, topography, topographic shading, permafrost occurrence, katabatic winds, monthly frequency of snow-free days) and vertical (depth and tree cover) processes driving microclimate variability⁴¹. Empirically-estimated coefficients were then used to assess temporal variation of microclimatic conditions between 2001–2005 and 2016–2020. The comparison of estimates from these two periods allowed measuring long-term annual and seasonal microclimate variations, provided estimates of the global-scale buffering effects of microclimate, and revealed that recently deglaciated habitats are globally experiencing a much stronger microclimatic change compared to other high-elevation environments.

Temperature modelling

The 175 microclimatic stations provided 706,810 temperature records in the period 2011–2021 (see Methods). Soil temperature was positively related to macroclimatic temperature, potential incoming solar radiation, frequency of snow-free days and distance from the glacier forefront (indicating a weak effect of katabatic winds), while it was negatively related to monthly cloud cover, occurrence of permafrost and tree cover (Fig. 2e-f). Statistically significant interactions showed that the increase of soil temperature at growing macroclimatic temperature was faster as the frequency of snow-free days increased (e.g. moving from spring to summer; Fig. 2a), while the interactions between snow-free days and both incoming radiation (Fig. 2b) and distance from the glacier were weaker. Higher cloud cover was associated with colder soil, with slightly stronger effects of cloud cover when the frequency of snow-free days was small (i.e. at the beginning of the season; Fig. 2c). Cloud cover also reduced the increase of soil temperature due to incoming solar radiation (Fig. 2d). Recorded temperature was also affected by the depth of burying, with loggers placed close to the surface recording higher temperatures (Fig. 2g). We did not detect an effect of the interaction between the monthly frequency of snow-free days and the distance from the glacier (Extended Data Table 1). The model provided a very good fit to the observations and explained a very high portion of microclimatic variations (R²_m = 0.75, R²_c = 0.87).

Downscaled macroclimate was the strongest driver of soil temperature, considering either variable importance scores (i.e. joint contribution of both additive and interactive terms; Extended Data Fig. 1a), or semi-partial R² (Extended Data Fig. 1b). Among the remaining predictors, solar radiation, monthly frequency of snow-free days, cloud cover and permafrost occurrence all explained a substantial portion of the total variance in soil temperature, either alone or interacting with each other, while the contribution of depth of burying, tree cover and distance from the glacier was small (Extended Data Fig. 1).

Extended Data Fig. 1

Variable contribution to the full model in terms of a) variable importance score (single predictors, measuring the joint contribution to both additive and interactive terms) and b) semi-partial R² (single terms). Error bars represent the 95% confidence intervals for the average estimate, obtained with 1,000 randomizations for each predictor (a) or 1,000 bootstrap replicates (b).

The evaluation of model transferability using a leave-one-out approach confirmed the robustness of our results. When any one of the glaciers was excluded, marginal and conditional R² values remained stable (average R²_m = 0.75, 2.5%-97.5% quantiles = 0.73–0.78; average R²_c = 0.87, 2.5%-97.5% quantiles = 0.87–0.89). Soil temperature values predicted during the cross-validation process were in good agreement with the observed ones (Fig. 2h), using both the unweighted set (coefficient of determination - R² = 0.73; mean absolute error - MAE = 1.33) and in the set adjusted to account for between-glacier differences in sample size and within-glacier differences in the frequency of snow-free days (weighted coefficient of determination - wR² = 0.82; weighted mean absolute error - wMAE = 1.44°C). When we predicted soil temperature using model coefficients averaged across the 26 leave-one-out models, we obtained nearly identical performance to the complete model (wR² = 0.85; wMAE = 1.30°C; Extended Data Table 2), confirming the predictive efficiency of our model.

To understand how adding predictors besides macroclimate improves the prediction of soil temperature, and the effect of downscaling macroclimate, we compared observed temperatures with i) those predicted by the complete model; ii) the downscaled macroclimate only; and the time-series of two widely used climate products: iii) TerraClimate⁴² and iv) Chelsa⁴³ (Extended Data Fig. 2). Our model outperformed both the traditional climate products and the downscaled macroclimate in predicting soil temperature, in terms of both variance explained (wR² = 0.83 vs 0.52–0.68) and mean absolute error (wMAE = 1.45°C vs 1.98 to 2.52°C; 1:1 line).

Extended Data Fig. 2: Comparison between the performance of our model and alternative approaches to the estimation of local temperature. Due to reduced temporal extent of the Chelsa dataset, the observations from 2020 were excluded from all the comparisons, and the recorded soil temperature was regressed against a) the predictions obtained from the leave-one-out approach, b) the downscaled macroclimate and the estimates obtained with the widely used climate products c) TerraClimate and d) Chelsa. The dashed black lines mark the perfect fit (1:1 line), while the red lines represent the fits from the weighted linear regression; shaded red areas represent the 95% confidence interval of the average estimates. To evaluate performances, the weighted coefficient of determination (wR²) and mean absolute error (wMAE) are provided.

Global projections

Building upon the high transferability of our modelling approach, we upscaled the model at the global scale for the periods 2001–2005 and 2016–2020. During this period, we detected substantial temperature increases in North America, the Andes and the higher latitudes of the Eastern Palaearctic, as well as in the European Alps and some areas of the Himalayas (Fig. 3a). When looking at different latitudinal bands (Fig. 3b), the pattern of temperature change showed clear differences between the Inter-tropical zone (23.44° S to 23.44° N), the Northern (23.44° to 72° N) and the Southern (23.44° to 60° S) Hemispheres. Temperature increase was particularly marked in the Inter-tropical zone and the Southern Hemisphere (mean ± sd: 0.68 ± 0.29 and 0.61 ± 0.45°C, respectively), with a generally higher increase in areas closer to glaciers (mean increase within 100 m from the glacier outline; Inter-tropical: 0.72 ± 0.33, Southern Hemisphere: 0.84 ± 0.62°C; increase 3 km from glaciers: 0.57 ± 0.26 and 0.59 ± 0.37°C, respectively; Fig. 3b). In the Northern Hemisphere, the change was smaller (0.41 ± 0.51°C), still temperature increase remained higher nearby glaciers compared to areas located 3 km away from the glacier (0.49 ± 0.70°C vs. 0.38 ± 0.40°C; Fig. 3b).

The analysis of seasonal trends provided comparable results, in terms of spatial patterns of temperature changes (Extended Data Fig. 3). In the mountain ranges of the Northern Hemisphere, the strongest temperature increase occurred during March-May and September-November, with a particularly intense increase in North America (Extended Data Fig. 3b-d). In the Southern Hemisphere, temperature changes were especially strong from September to May (Extended Data Fig. 3). Conversely, in the Inter-tropical zone temperature change was homogeneously distributed throughout the year, owing to the reduced effect of seasonality. A decrease of changes with increasing distance from the glacier was evident during the period September-February (Extended Data Fig. 4a and d). Seasonal trends confirmed a stronger temperature increase nearby glaciers (Extended Data Fig. 4); this effect was evident for all latitudinal bands, and more evident during warm seasons.

Extended Data Fig. 3

Seasonal trends of temperature change between 2016–2020 and 2001–2005. Per-cell average changes in soil temperature during a) December-February, b) March-May, c) June-August and d) September-November; the dashed horizontal lines identify the Tropics, while the continuous line indicates the Equator.

We estimated the potential for microclimatic buffering as the ratio between the current spatial variability and the temporal change experienced during the last 20 years. The sign of this relationship returns the direction of the temporal change (i.e. either temperature increase or decrease), while its absolute value measures the buffering potential (e.g. values of + 2 or -2 indicate spatial variability twice the temporal variation, given increasing or decreasing temperatures, respectively). In almost all cases, the fine-scale spatial variability of soil temperature within 250 m was larger than the temporal change, suggesting that it can play a relevant role for microclimatic buffering (Fig. 3c). The majority of buffering values (57.56–61.13%, depending on the latitudinal band and the distance from the glacier) had values > 1 and ≤ 10, indicating that in the last 20 years spatial variability was one-to-ten times larger than the temporal temperature change. When considering all values ≤ -1 or > 1 (i.e. looking at all the sites potentially guaranteeing buffering), percentages ranged between 89.42 to 98.04%.

The duration of the snow-free season estimated from satellite images increased between 2001–2005 and 2016–2020. At the global scale, the mean increase of season duration was 9.44 days (sd: 21.48 days; Fig. 3d), but the increase was larger nearby glaciers (mean ± sd: 16.42 ± 25.35 days) compared to areas located at 3 km from a glacier (mean ± sd: 5.28 ± 18.09 days). In the Inter-tropical zone, the effect of distance from glaciers was particularly marked. Here, almost no change in season duration occurred in sites located more than 1 km away from a glacier, probably because in the tropics areas far from glaciers are almost constantly without snow.

Extended Data Fig. 4

Seasonal trends of temperature change between the periods 2016–2020 and 2001–2005. For each latitudinal band and distance class, violin plots summarize temperature changes in soil temperature during a) December-February, b) March-May, c) June-August and d) September-November. Yellow dots mark the median value for each series, while yellow lines the first and third quartiles.

Understanding the effects of climate change on high-elevation ecosystems is pivotal to predict the future of these threatened environments. Accurate microclimatic information is crucial to identify the conditions that are effectively experienced by living organisms, as changes in microclimate strongly influence local distribution and survival of individuals²⁴. At the same time, the extreme environmental heterogeneity of mountain habitats, mainly generated by the altitudinal gradients, the complex topography and the variability of vegetation (from mature forests to grasslands, peatlands and tundra to bare soils) determines patchy microhabitats with a wide range of microclimatic conditions in relatively small areas, potentially buffering the severity of warming impacts on populations^35,36,37.

As expected, temporal changes in microclimate are tightly linked to climate trends at the regional or global scale, with macroclimate playing the major role in driving local temperatures. Our results highlighted a generalized increase in soil temperature between 2001 and 2020 at all latitudes and distances from the glacier front, with the presence of clear seasonal trends. In both hemispheres, microclimate variation was stronger (> 0.5°C) during warm seasons, while it was reduced during cold months. In the same period, we recorded an increase in the annual duration of the snow-free season, confirming at the global scale the results of regional analyses³¹. Both these changes were particularly evident nearby glaciers, where the shrinkage of ice with the consequent reduction of its cooling effect amplifies the temperature rise. The increasing number of days with snow-free terrains deeply influences soil temperature, mainly through the interaction with other drivers of temperature (Fig. 2). The absence of snow cover increases heat exchanges between air and soil, and the lower albedo facilitates absorption of solar radiation, resulting in steeper relationships²⁵. This amplifies the temperature increase along elevational gradients⁴⁴, and likely has major impacts on the whole ecosystem, such as an increase in vegetation productivity, and a change of biotic communities^31,33.

Owing to the mechanism of elevation-dependent warming, temperature increase is faster in mountain areas than in surrounding lowlands^44,45. Such accelerated warming poses strong challenges to the mountain ecosystems, potentially leading to species altitudinal migrations, phenological changes and mismatches between different components of the ecosystem⁴⁶. However, the impact of increasing soil temperatures and duration of snow-free season on local alpine biota may be partially counterbalanced by the spatial variability of microclimate conditions. For example, Maclean⁴⁷ recorded differences in temperature of almost 20°C across a four hectares study area, corresponding to differences in temperature across entire continents. When analysing the potential for microclimate buffering, we found conspicuous variations in soil temperature, with a mean value of about 2°C within a limited spatial range (250 m), i.e. one-to-ten times the recorded change during the last 20 years. Such spatial variability in microclimate conditions has key effects on local communities^33,48, and might allow individuals and even communities to withstand, at least temporarily, the effects of climate warming by modifying their distribution over relatively limited distances. Nevertheless, fine-scale heterogeneity is probably not enough to buffer warming patterns expected to occur in the long term, as several climate change scenarios suggest that most of mountain regions will experience a warming > 4°C by the end of the century².

Differences in soil moisture may influence local temperature, with moister soils having higher thermal inertia, and soil water content buffering temperature variations⁴⁹. Differences in water availability also influence the local distribution and survival of individuals⁵⁰, nitrogen mineralization⁵¹, the fluxes of carbon across the soil-vegetation-atmosphere interface²¹ and the development of alpine ecosystems^12,19. The combined effects of differences in soil moisture and temperature may result in microclimate patterns even more complex than those generated by temperature alone; still, the lack of high-resolution data on soil properties (e.g. texture, composition) hampers the modelling of soil water content at regional and global scale⁴⁰. As far as we know, the dataset we assembled is the most complete collection of soil temperature recordings in high-mountain areas. Despite we tried to cover areas at different latitudes and with diverse climatic conditions, we acknowledge that data used to develop our model are not truly global, thus may not be fully representative of the conditions in other mountain ranges with very different climatological characteristics. For instance, in some regions glaciers currently are not retreating, owing to the combination of locally stationary temperatures, increasing precipitation and / or heavy debris cover (e.g. Karakoram⁵²); in others they are considerably retreating, but following local dynamics (e.g. the southern Himalayas, conditioned by the Indian monsoon). A further extension of the dataset and its implementation with other initiatives (e.g. 41) might improve the global representation of the overall high-elevation soil temperature dynamics, still our analysis provide much better information on the microclimate of high mountain environments than any currently available product (Extended Data Fig. 2), and enables an unprecedented view of the fine-scale heterogeneity of climate change.

During the last decades we boosted our understanding of the drivers of microclimate²⁴. A huge amount of information is already available at the macro- and meso-scales for long time-series (e.g. monthly soil moisture, shortwave radiation or temperatures from TerraClimate), and other can be retrieved using remotely-sensed data. However, the effective and consistent modelling of microclimatic conditions on large areas still requires important efforts. Blending process-based models (e.g. 53) with data-driven empirical models (e.g. 21) and assimilating the data flow produced by remote sensing could be a first step to the construction of a “digital twin” of alpine ecosystems, which will be fundamental to model and understand relationships between mountain species and their environment, and to quantify their responses to climate change.

Acknowledgements

S.M., M.C., M.G., R.A., R.S.A., F.G., F.P., W.T. and G.F.F. are funded by the European Research Council under the European Community’s Horizon 2020 Programme, grant agreement no. 772284 (‘IceCommunities-Reconstructing community dynamics and ecosystem functioning after glacial retreat’).

Author Contributions

S.M. and G.F.F. conceived and designed the study. S.M., A.Z, M.C, M.G., R.A., R.S.A., F.G, F.P. and G.F.F. collected and validated field data. S.M. wrote the codes, conducted the statistical analyses and produced the illustrations with major inputs from G.F.F., A.P and W.T. All authors contributed to writing the article.

Data availability

Soil temperature data will be available on FigShare upon article acceptance. Codes are provided in Supplementary Software 1.

Competing interests

The authors declare no competing interests.

Supplementary Information

Supplementary Table 1 and Supplementary Software 1

Nogués-Bravo, D., Araújo, M. B., Errea, M. P. & Martínez-Rica, J. P. Exposure of global mountain systems to climate warming during the 21st Century. Glob. Environ. Change 17, 420–428 (2007).
Hock, R. et al. High Mountain Areas. In: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate (Intergovernmental Panel on Climate Change, 2019).
Zimmer, A., Beach, T., Klein, J. A. & Recharte Bullard, J. The need for stewardship of lands exposed by deglaciation from climate change. WIREs Clim. Change 13, e753 (2022).
Rahbek, C. et al. Building mountain biodiversity: Geological and evolutionary processes. Science 365, 1114–1119 (2019).
Rahbek, C. et al. Humboldt’s enigma: What causes global patterns of mountain biodiversity?. Science 365, 1108–1113 (2019).
Körner, C. et al. A global inventory of mountains for bio-geographical applications. Alp. Bot. 127, 1–15 (2017).
Hoorn, C., Perrigo, A. & Antonelli, A. Mountains, climate and biodiversity (John Wiley & Sons, 2018).
Hugonnet, R. et al. Accelerated global glacier mass loss in the early twenty-first century. Nature 592, 726–731 (2021).
Matthews, J.A. The ecology of recently-deglaciated terrain: a geoecological approach to glacier forelands (Cambridge University Press, 1992).
Beniston, M. et al. The European mountain cryosphere: a review of its current state, trends, and future challenges. The Cryosphere 12, 759–794 (2018).
Cauvy-Fraunié, S. & Dangles, O. A global synthesis of biodiversity responses to glacier retreat. Nat. Ecol. Evol. 3, 1675–1685 (2019).
Ficetola, G. F. et al. Dynamics of ecological communities following current retreat of glaciers. Annu. Rev. Ecol. Evol. Syst. 52, 405–426 (2021).
Chapin, F. S., Walker, L. R., Fastie, C. L. & Sharman, L. C. Mechanisms of primary succession following deglaciation at Glacier Bay, Alaska. Ecol. Monogr. 64, 149–175 (1994).
Khedim, N. et al. Topsoil organic matter build-up in glacier forelands around the world. Glob. Change Biol. 27, 1662–1677 (2021).
Erschbamer, B., Niederfriniger Schlag, R. & Winkler, E. Colonization processes on a central Alpine glacier foreland. J. Veg. Sci. 19, 855–862 (2008).
Jenny, H. Factors of soil formation: a system of quantitative pedology (Courier Corporation, 1994).
Walker, L. R., Wardle, D. A., Bardgett, R. D. & Clarkson, B. D. The use of chronosequences in studies of ecological succession and soil development. J. Ecol. 98, 725–736 (2010).
Cannone, N., Diolaiuti, G., Guglielmin, M. & Smiraglia, C. Accelerating climate change impacts on alpine glacier forefield ecosystems in the European Alps. Ecol. Appl. 18, 637–648 (2008).
Vater, A. E. & Matthews, J. A. Succession of pitfall-trapped insects and arachnids on eight Norwegian glacier forelands along an altitudinal gradient: Patterns and models. The Holocene 25, 108–129 (2015).
Fickert, T., Grüninger, F. & Damm, B. Klebelsberg revisited: did primary succession of plants in glacier forelands a century ago differ from today?. Alp. Bot. 127, 17–29 (2017).
Magnani, M. et al. Drivers of carbon fluxes in Alpine tundra: a comparison of three empirical model approaches. Sci. Total Environ. 732, 139139 (2020).
Magnani, M., Baneschi, I., Giamberini, M., Raco, B. & Provenzale, A. Microscale drivers of summer CO2 fluxes in the Svalbard High Arctic tundra. Sci. Rep. 12, 763 (2022).
Dobrowski, S. Z. A climatic basis for microrefugia: the influence of terrain on climate. Glob. Change Biol. 17, 1022–1035 (2011).
Bramer, I. et al. Advances in monitoring and modelling climate at ecologically relevant scales. Adv. Ecol. Res. 58, 101–161 (2018).
Niittynen, P. & Luoto, M. The importance of snow in species distribution models of arctic vegetation. Ecography 41, 1024–1037 (2018).
Lembrechts, J. J. et al. Global maps of soil temperature. Glob. Change Biol. 28, 3110–3144 (2022).
Thornthwaite, C. W. A charter for climatology. WMO Bulletin 2, 40–46 (1953).
Pepin, N. C. & Lundquist, J. D. Temperature trends at high elevations: patterns across the globe. Geophys. Res. Lett. 35, L14701 (2008).
Hågvar, S. et al. Ecosystem birth near melting glaciers: a review on the pioneer role of ground-dwelling arthropods. Insects 11, 644 (2020).
Callaghan, T. V. et al. The changing face of Arctic snow cover: A synthesis of observed and projected changes. Ambio 40, 17–31 (2011).
Rumpf, S. B. et al. From white to green: Snow cover loss and increased vegetation productivity in the European Alps. Science 376, 1119–1122 (2022).
Niittynen, P., Heikkinen, R. K. & Luoto, M. Decreasing snow cover alters functional composition and diversity of Arctic tundra. Proc. Natl Acad. Sci. USA 117, 21480–21487 (2020).
Körner, C. Alpine plant life: functional plant ecology of high mountain ecosystems (Springer Nature, 2021).
De Frenne, P. et al. Forest microclimates and climate change: Importance, drivers and future research agenda. Glob. Change Biol. 27, 2279–2297 (2021).
Lenoir, J., Hattab, T. & Pierre, G. Climatic microrefugia under anthropogenic climate change: implications for species redistribution. Ecography 40, 253–266 (2017).
Suggitt, A. J. et al. Extinction risk from climate change is reduced by microclimatic buffering. Nat. Clim. Change 8, 713–717 (2018).
De Frenne, P. et al. Global buffering of temperatures under forest canopies. Nat. Ecol. Evol. 3, 744–749 (2019).
Feldmeier, S. et al. Shifting aspect or elevation? The climate change response of ectotherms in a complex mountain topography. Diver. Distrib. 26, 1483–1495 (2020).
Maclean, I. M., Suggitt, A. J., Wilson, R. J., Duffy, J. P. & Bennie, J. J. Fine-scale climate change: modelling spatial variation in biologically meaningful rates of warming. Glob. Change Biol. 23, 256–268 (2017).
Lembrechts, J. J., Nijs, I. & Lenoir, J. Incorporating microclimate into species distribution models. Ecography 42, 1267–1279 (2019).
Lembrechts, J. J. et al. SoilTemp: A global database of near-surface temperature. Glob. Change Biol. 26, 6616–6629 (2020).
Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Scientific Data 5, 170191 (2018).
Karger, D. N. et al. Climatologies at high resolution for the earth’s land surface areas. Scientific Data 4, 170122 (2017).
Pepin, N. et al. Elevation-dependent warming in mountain regions of the world. Nat. Clim. Change 5, 424–430 (2015).
Pepin, N. C. et al. Climate changes and their elevational patterns in the mountains of the world. Rev. Geophys. 60, 2020RG000730 (2022).
Vitasse, Y. et al. Phenological and elevational shifts of plants, animals and fungi under climate change in the European Alps. Biol. Rev. 96, 1816–1835 (2021).
Maclean, I. M. Predicting future climate at high spatial and temporal resolution. Glob. Change Biol. 26, 1003–1011 (2020).
Viterbi, R., Cerrato, C., Bionda, R. & Provenzale, A. Effects of temperature rise on multi-taxa distributions in mountain ecosystems. Diversity 12, 210 (2020).
Ashcroft, M. B. & Gollan, J. R. Moisture, thermal inertia, and the spatial distributions of near-surface soil and air temperatures: understanding factors that promote microrefugia. Agric.For. Meteorol. 176, 77–89 (2013).
Marta, S., Brunetti, M., Manenti, R., Provenzale, A. & Ficetola, G. F. Climate and land-use changes drive biodiversity turnover in arthropod assemblages over 150 years. Nat. Ecol.Evol. 5, 1291–1300 (2021).
Raffl, C., Mallaun, M., Mayer, R. & Erschbamer, B. Vegetation succession pattern and diversity changes in a glacier valley, Central Alps, Austria. Arct. Antarct. Alp. Res. 38, 421–428 (2006).
Farinotti, D., Immerzeel, W. W., de Kok, R. J., Quincey, D. J. & Dehecq, A. Manifestations and mechanisms of the Karakoram glacier Anomaly. Nat. Geosci. 13, 8–16 (2020).
Baudena, M., d'Andrea, F. & Provenzale, A. A model for soil-vegetation‐atmosphere interactions in water‐limited ecosystems. Water Resour. Res. 44, 2008WR007172 (2008).
Raup, B. et al. The GLIMS geospatial glacier database: a new tool for studying glacier change. Glob. Planet. Change 56, 101–110 (2007).

Temperature data and predictors

Data on soil temperature were collected at 175 sites from 26 glacier forelands located in the Svalbard archipelago (Norway; 2 forelands), European Alps (Italy, Austria, Switzerland and France; 21 forelands), and the Andes (Peru; 3 forelands), between 15^th July 2011 and 24^th August 2021 (Figure 1). In each foreland, 1 to 16 devices (mean: 6.7; sd: 4.3; total: 175) were buried with no shielding at 5, 10 or 15 cm (mean: 9.2; sd: 3.0) below soil surface. The distance between devices within the same foreland ranged from 0.1 to 1798 m (mean: 149; sd: 223). Devices recorded temperatures for 33 to 763 days (mean: 501.9; sd: 191.7) with a recording frequency varying between 6 and 480 recordings/day (mean: 11.1; sd: 36.5). The total dataset was composed by 706,810 recordings; device model, recording parameters and burying depth showed some difference depending on the foreland (see Supplementary Table 1 for technical specifications). Dataset cleansing involved the removal of i) all data collected after 31^st December 2020, and ii) months sampled for less than 90% of time. Monthly averages were calculated based on the remaining data. The final dataset was composed of 2,203 monthly average temperatures from 26 glacier forelands, ranging between August 2011 and December 2020. For each glacier foreland, the region of interest (ROI) was defined as the extent enclosing all the sampling stations, with a 750 m buffer as this enabled to include areas from the glacier tongue to downstream areas; a larger buffer (1500 m) was set up for Morteratsch and Dammagletscher forelands in order to include part of the glacier tongue within the ROI.

Macroclimate information for the sampled months and years was retrieved from the medium-resolution climate product TerraClimate⁴² (150 arcsec). Monthly mean temperatures were calculated as the midpoint between TerraClimate monthly minimum and maximum temperatures. To account for the adiabatic decrease of temperature we applied the standard and fixed environmental lapse rate of -0.0065 °C/m⁵⁵. The high-resolution digital elevation data needed for downscaling macroclimate were retrieved from the ASTER GDEM v3 (resolution: 1 arcsec, i.e. approx. 30 m at the equator; latitudinal extent: 82° N to 83° S) via the NASA Earthdata interface (https://doi.org/10.5067/ASTER/ASTGTM.003; last accession on 24^th March 2020). The high-resolution temperature surfaces obtained downscaling TerraClimate simply represented the fine-scale variability of macroclimate related to altitudinal differences, rather than the effective topo- or micro-climates.

To account for the effect of topography, we calculated the incoming solar radiation using the solarindex function³⁹, with slope and aspect data retrieved from ASTER GDEM v3. Given a location and time, solarindex returns the proportion of direct beam radiation intercepted by a surface, taking into account solar altitude and azimuth, and topographic shading. For any given month and year, we i) calculated the index iteratively at hourly intervals (0 to 23) for three days per month (the 5^th, 15^th and 25^th)³⁶; ii) summed the hourly estimates of each day, to obtain the daily cumulative incoming radiation; and iii) averaged daily estimates to obtain the average cumulative daily radiation for the month of interest. Cloud cover reduces the incoming solar radiation intercepted by surfaces. Coarse-grained (30 arcsec) MODIS-derived monthly frequency of clouds, averaged over the period 2000-2014⁵⁶ was retrieved from EarthEnv (http://www.earthenv.org/cloud; accessed on 27^th September 2021) and bilinearly interpolated to generate cloudiness profiles at 1 arcsec resolution.

Snow cover (i.e. the presence of snow on the ground) causes local soil temperatures strongly decoupled from regional climate²⁶, due to the insulating effect of the snowpack. To quantify this thermal effect on soil temperatures, for each month we assessed the proportion of days in which the device was under the snow. When a device is under the snow, it shows a very limited daily variation. We tested several values of diurnal range, assuming that a sensor was under snow when it showed a range below different threshold values (0.5, 1, 1.5, 2, 2.5 and 3°C), and calculated the corresponding number of days with no snow on the ground for each of these thresholds. These estimates were compared with the ones obtained, for the same months and years, using the NDSI-derived coarse-grained MODIS Terra 500m daily fractional snowcover⁵⁷ (available at https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD10A1). Fractional snowcover was converted to snow occurrence using the conservative threshold of 40%⁵⁸, and the monthly frequency of snow-free days was estimated. Whatever the threshold, we found an excellent match between the number of days under the snow estimated from sensor and from the MODIS data (Pearson’s r ranging from 0.918 to 0.942); to calculate the proportion of snow-free days, the highest correlation was obtained for a threshold value of 1.5°C (r: 0.942). Therefore the proportion of snow-free days was calculated on the basis of this threshold.

The presence of nearby glaciers represents a further driver of local temperature, for instance because of katabatic winds. To account for this cooling effect, we measured the distance between each sampling station and the glacier front, under the assumption of decreasing cooling effects with distance⁵⁹. For each glacier, the most recent outline was retrieved from Marta et al.⁶⁰, checked against glacier position between 2015 and 2019 using USGS Landsat 8 imagery (available at https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_TOA) and eventually updated to include other nearby glaciers or more recent outlines. Outlines were transformed to polygons, rasterized, and distance maps at 30 m resolution were calculated using the function gridDistance from the raster R package⁶¹. Three categorical variables were calculated to account for the effect of i) tree cover, ii) permafrost occurrence and iii) differences in the depth of logger burying. Tree shading can reduce soil temperature, thus we used the 30-m resolution Hansen Global Forest Change v1.8⁶² (available at https://developers.google.com/earth-engine/datasets/catalog/UMD_hansen_global_forest_change_2020_v1_8) to assess if a logger was in a tree-covered area or not. Data on permafrost extent were retrieved from Gruber⁶³; given the coarse-grained resolution of this dataset (30 arcsec), we calculated the spatially weighted mean over each ROI after disaggregating at the 30 m resolution. To convert mean permafrost extent over the ROI to occurrence of continuous / extensive discontinuous permafrost we used the conservative 0.9 threshold; i.e. we considered permafrost present over the ROI when the probability of occurrence was ≥ 0.9.

Model calibration and testing

Monthly-averaged observed soil temperature (soilT) was modelled using linear mixed models (LMM). As independent variables we used macroclimate (mT), daily cumulative potential incoming solar radiation (rad), monthly percent cloud cover (clo), monthly frequency of snow-free days (sfd), distance from glacier forefront (dg), tree cover (tc), permafrost occurrence (pf) and depth of burying (d). To account for the effects of a varying frequency of snow-free days (sfd) on mT, rad, clo and dg, as well as for the effect of clo on rad, interactive terms were added. To include geographical factors not explicitly accounted for by the selected set of predictors, we additionally incorporated a random intercept on glacier (1|gl). Consequently, the full model takes the form:

soilT ~ mT + rad + clo + sfd + dg + additive effects

sfd:mT + sfd:rad + sfd:clo + sfd:dg + sfd - interactive effects

rad:clo + clo - interactive effects

tc + pf + d + intercept corrections

1|gl random intercept

Winter snowpack decouples air and soil temperatures causing no relationship between soilT and several predictors (e.g. air temperature, solar radiation, cloud cover) during seasons with snow. To remove the effects of this decoupling, while reconstructing soilT during the snow-free season, we i) classified sfd in 10 intervals between 0 and 100%; ii) run iteratively the model retaining only records with sfd > 10%, 20%,..., and iii) plotted fitted values vs residuals to evaluate the residual structure at each step. With sfd > 30% the effect of decoupling was almost completely erased. Consequently, all the months with sfd ≤ 30% were discarded, and the resulting dataset included 1,258 monthly average temperatures from 26 glacier forelands. Before running the final model, dg was square-root transformed to linearize the relationship with the response variable and all the continuous predictors were scaled to zero mean and unit variance. Linear mixed models were run using the lme4⁶⁴ and lmerTest⁶⁵ R packages. Model residuals approximated a normal distribution (Shapiro-Wilk test; W = 0.996), and the variance inflation factor was low (GVIFadj_max= 1.81, including interaction terms), indicating that multicollinearity did not pose major issues. Model performances were evaluated using Nakagawa and Schielzeth⁶⁶ R², as implemented in MuMIn R package⁶⁷. The amount of variance explained by single model terms was quantified by calculating the semi-partial R² using the partR2 R package⁶⁸, with 1,000 bootstrap replicates. partR2 iteratively removes predictors, and compares the change in variance of the linear predictor to the variance explained by the full model; higher the difference between the two values, higher the amount of variance explained uniquely by a given predictor. We followed Thuiller et al.⁶⁹ to account for the overall effect of single predictors (i.e. considering their joint contribution to both additive and interactive terms). Each predictor was randomized 1,000 times, and the predictions obtained using original and randomized datasets were compared via the Pearson’s correlation coefficient (r). Strong correlations indicate that randomizations had little effect on model performances; for each permutation, variable importance was finally expressed as 1-r.

To evaluate model transferability, we followed a leave-one-out approach. We iteratively run the full model, retaining all glaciers except one as the training set; the estimated fixed coefficients were then used to predict the expected temperature for the excluded glacier. Differences in sample size between glaciers, as well as those in the monthly frequency of snow-free days within each glacier may inflate agreement scores. To account for these differences, we thus conservatively downweighted each observation, so that observations of each glacier received a weight = 1/M, where M is the number of recordings in any given sfd category available for that glacier. The agreement between observed and predicted temperatures was measured using i) the coefficient of determination form a weighted linear regression (wR²), and ii) the weighted mean absolute error (wMAE).

To confirm that temperatures estimated by model approximate the actual temperature better than other already available products, we also compared observed temperatures to: the ones predicted by the model, the downscaled macroclimate (mT), the time-series of TerraClimate⁴² (resolution: 150 arcsec) and Chelsa⁴³ (resolution: 30 arcsec). For each observation, we extracted climate data for the corresponding year and month, after excluding the observations from 2020 (Chelsa dataset being limited to 2019), and calculated wR² and wMAE.

Global projection of soil temperature

Obtaining high-resolution estimates of soil temperature in glacier forelands, at the global scale and in several periods, allows estimating soil microclimate variability and temporal variation, measuring the impacts of climate change on microclimate and the potential for microclimate buffering. The aim of this analysis was to assess the variation of microclimate during the last decades, thus we compared microclimate between the periods 2001-2005 and 2016-2020. We used the mean coefficients obtained from the leave-one-out analysis to generate predictions of soil temperature at the global scale, using Google Earth Engine and the rgee R package⁷⁰. Due to data availability, the analysis was spatially constrained between 60° S and 72° N. We focused on proglacial landscapes, thus we limited projections to within 3 km from glacier outlines. It is worth noting that some differences exist between the calibration and global projection for: the digital elevation products, the approach to glacier outline identification, the definition of the monthly duration of the monthly frequency of snow-free days and the calculation of potential incoming solar radiation.

Digital elevation data are needed for downscaling macroclimate and for calculating the daily cumulative potential incoming solar radiation. For the global projection, we used a coarser resolution (90 m instead of 30 m) to limit computation time. From 60° S to 60° N we used the 90 m resolution composite of Shuttle Radar Topographic Mission v4⁷¹ (available at https://developers.google.com/earth-engine/datasets/catalog/CGIAR_SRTM90_V4), while from 60° to 72° N we used the Global Multi-resolution Terrain Elevation Data 2010⁷² (available at https://developers.google.com/earth-engine/datasets/catalog/USGS_GMTED2010), given the SRTM model was not available above 60 °N. Monthly mean temperature was calculated from TerraClimate (https://developers.google.com/earth-engine/datasets/catalog/IDAHO_EPSCOR_TERRACLIMATE). Monthly mean temperatures were obtained by averaging monthly minimum and maximum values across each five-year period, and downscaled to 90 m resolution following the same approach used for the calibration data. To obtain information on the potential incoming solar radiation, the original solarindex function was re-coded to be launched directly in GEE via the rgee interface (see Supplementary Software 1). For each cell, daily cumulative incoming radiation was estimated for the 15^th day of each month in the years 2003 (for 2001-2005) and 2018 for (2016-2020), and considered representative of the whole month. Monthly frequency of cloud cover was uploaded in GEE, and bilinearly-interpolated at the 90 m resolution. The monthly frequency of snow-free days was calculated using the NDSI-derived daily fractional snowcover as detailed in the previous section. For each month, we estimated the percent snow occurrence using monthly values averaged over each five-year period and bilinearly-interpolated at the 90 m resolution. All cells with sfd values ≤ 30% were masked (i.e. excluded from the analysis). To account for distance from the glacier, we used glacier outlines of the GLIMS dataset⁵⁴ (available at https://developers.google.com/earth-engine/datasets/catalog/GLIMS_current). Glaciers may have been retreating between 2001-2005 and 2016-2020; consequently, for each period and glacier (“glac_id” field), we selected the outline with the temporally closer source image (“src_date” field), and calculated distances according to those positions. Permafrost extent was uploaded in GEE, and bilinearly-interpolated at the 90 m resolution, while tree cover was aggregated at the same 90 m resolution. In projections, we estimated soil temperature at 5 cm depth (d = 5). Despite the methodological differences, incoming radiation and temperatures estimated with the global model (90 m resolution) showed excellent agreement with the ones obtained by the 30 m resolution (Pearson’s r = 0.942 and 0.924, respectively).

Maps of predicted soil temperatures at 2001-2005 and 2016-2020 pose some problems in handling and obtaining summary statistics at the global scale (2 periods × 12 months × 6.628 × 10¹⁰ pixels; approximate size ≈ 4.7 TB). To overcome these limitations and obtain a spatially unbiased representation of microclimate variability and variation, instead than using all the cells we subsampled them using a stratified grid sampling by i) building a regular 50 × 50 km grid (Mollweide projection; ESRI:54009); ii) retaining all the grid cells containing glaciers or within 3 km from glacier outlines (2,604 cells), and iii) defining five classes of distance from the most recent glacier outline (0-100, 400-600, 900-1,100, 1,900-2,100 and 2,900-3,100 m). The most recent glacier outline was the one used for calculating the distances for the 2016-2020 projection. Within each cell, we randomly sampled 10 points, two for each distance class. The resulting dataset was composed of 26,040 points, each associated with 12 × 2 (2001-2005 and 2016-2020) measures of monthly soil temperature. After removing points with missing temperature estimates in one or both periods, and cells with < 9 points, the final dataset was composed of 22,415 records from 2,274 cells. For this set of points, we extracted the monthly average temperature and annual duration of the snow-free season for the two periods. Based on temperature data, we calculated both annual and seasonal (dec-feb; mar-may; jun-aug and sep-nov) microclimate variation (DT) between the two periods (DT = T_2016-2020- T_2001-2005).

Short-distance movement of individuals might allow buffering the severity of warming impacts on populations, if suitable climatic conditions occur nearby³⁶. To understand the potential for microclimate buffering of proglacial environments, we compared the recorded microclimate variation between 2001-2005 and 2016-2020 (DT) to the spatial variability of soil temperatures. The spatial variability of microclimate was calculated as the 80% inter-percentile range within a 250 m buffer (T_var). Due to computing limitations, the analysis was restricted to five points per cell, one for each class of distances from the glacier, considering the average microclimate (mean annual temperature) of 2016-2020. Microclimate buffering potential (T_bp) was calculated as: T_bp = T_var/DT. This formula allows measuring both the direction of the change and the buffering potential, as it retains the sign from DT (e.g. positive values indicate temperature increase), but returns (absolute) values ≥ 1 (|T_bp| ≥ 1) when the spatial microclimate variability is larger than the temporal microclimate variation.

Methods-only References

55. Barry, R. G. Mountain Weather and Climate (Cambridge University Press, 2008).

56. Wilson, A. M. & Jetz, W. Remotely sensed high-resolution global cloud dynamics for predicting ecosystem and biodiversity distributions. PLoS Biol. 14, e1002415 (2016).

57. Hall, D. K., Salomonson, V. V. & Riggs, G. A. MODIS/Terra Snow Cover Daily L3 Global 500m Grid - Version 6 (NASA National Snow and Ice Data Center - Distributed Active Archive Center, 2016).

58. Xiao, X. et al. Estimating fractional snow cover from passive microwave brightness temperature data using MODIS snow cover product over North America. The Cryosphere 15, 835-861 (2021).

59. Shaw, T. E. et al. Distributed summer air temperatures across mountain glaciers in the south-east Tibetan Plateau: temperature sensitivity and comparison with existing glacier datasets. The Cryosphere 15, 595-614 (2021).

60. Marta, S. et al. The retreat of mountain glaciers since the Little Ice Age: a spatially explicit database. Data 6, 107 (2021).

61. Hijmans, R. J. raster: Geographic Data Analysis and Modeling (R package version 3.4-5, 2020).

62. Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. Science 342, 850-853 (2013).

63. Gruber, S. Derivation and analysis of a high-resolution estimate of global permafrost zonation. The Cryosphere 6, 221-233 (2012).

64. Bates, D., Mächler, M., Bolker, B. & Walker S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1-48 (2015).

65. Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest Package: Tests in Linear Mixed Effects Models. J. Stat. Softw. 82, 1-26 (2017).

66. Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R² from generalized linear mixed‐effects models. Methods Ecol. Evol. 4, 133-142 (2013).

67. Barton, K. MuMIn: Multi-Model Inference (R package version 1.43.17, 2020).

68. Stoffel, M. A., Nakagawa, S. & Schielzeth, H. partR2: Partitioning R² in generalized linear mixed models. PeerJ 9, e11414 (2021).

69. Thuiller, W., Lafourcade, B., Engler, R. & Araújo, M. B. BIOMOD–a platform for ensemble forecasting of species distributions. Ecography 32, 369-373 (2009).

70. Aybar, C., Wu, Q., Bautista, L., Yali, R. & Barja, A. rgee: An R package for interacting with Google Earth Engine. J. Open Source Softw. 5, 2272 (2020).

71. Jarvis, A., Reuter, H. I., Nelson, A. & Guevara, E. Hole-filled SRTM for the globe - Version 4 (CGIAR-CSI SRTM 90m Database, 2008).

70. USGS. Global Multi-resolution Terrain Elevation Data (U.S. Geological Survey, 2010).

There is NO Competing Interest.

Martaetal.SupplementarySoftware1.txt
Martaetal.SupplementaryTable1.csv
ExtendedDataTables.docx
floatimage3.png
Extended Data Figure 1: Variable contribution to the full model in terms of a) variable importance score (single predictors, measuring the joint contribution to both additive and interactive terms) and b) semi-partial R² (single terms). Error bars represent the 95% confidence intervals for the average estimate, obtained with 1,000 randomizations for each predictor (a) or 1,000 bootstrap replicates (b).
floatimage4.png
Extended Data Figure 2: Comparison between the performance of our model and alternative approaches to the estimation of local temperature. Due to reduced temporal extent of the Chelsa dataset, the observations from 2020 were excluded from all the comparisons, and the recorded soil temperature was regressed against a) the predictions obtained from the leave-one-out approach, b) the downscaled macroclimate and the estimates obtained with the widely used climate products c) TerraClimate and d) Chelsa. The dashed black lines mark the perfect fit (1:1 line), while the red lines represent the fits from the weighted linear regression; shaded red areas represent the 95% confidence interval of the average estimates. To evaluate performances, the weighted coefficient of determination (wR²) and mean absolute error (wMAE) are provided.
Onlinefloatimage5.png
Onlinefloatimage6.png
Extended Data Figure 3: Seasonal trends of temperature change between 2016-2020 and 2001-2005. Per-cell average changes in soil temperature during a) December-February, b) March-May, c) June-August and d)September-November; the dashed horizontal lines identify the Tropics, while the continuous line indicates the Equator.
Onlinefloatimage7.png
Extended Data Figure 4: Seasonal trends of temperature change between the periods 2016-2020 and 2001-2005. For each latitudinal band and distance class, violin plots summarize temperature changes in soil temperature during a) December-February, b)March-May, c) June-August and d) September-November. Yellow dots mark the median value for each series, while yellow lines the first and third quartiles.

Download PDF

Journal Publication

published 31 Aug, 2023

Read the published version in Nature Communications →

Version 1

posted

You are reading this latest preprint version

Heterogeneous changes of soil microclimate in high mountains and glacier forelands

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Results

Temperature modelling

Global projections

Discussion

Declarations

References

Methods

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1