*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 15th July 2011 and 24th 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 31st 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 TerraClimate42 (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/m55. 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 24th 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* function39, 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 5th, 15th and 25th)36; 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-201456 was retrieved from EarthEnv (http://www.earthenv.org/cloud; accessed on 27th 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 climate26, 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 snowcover57 (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%58, 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 distance59. For each glacier, the most recent outline was retrieved from Marta et al.60, 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 package61. 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.862 (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 Gruber63; 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*64 and *lmerTest*65 R packages. Model residuals approximated a normal distribution (Shapiro-Wilk test; W = 0.996), and the variance inflation factor was low (GVIFadjmax = 1.81, including interaction terms), indicating that multicollinearity did not pose major issues. Model performances were evaluated using Nakagawa and Schielzeth66 *R*2, as implemented in *MuMIn* R package67. The amount of variance explained by single model terms was quantified by calculating the semi-partial *R*2 using the *partR2* R package68, 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.69 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 (w*R*2), 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 TerraClimate42 (resolution: 150 arcsec) and Chelsa43 (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 w*R*2 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 package70. 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 v471 (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 201072 (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 15th 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 dataset54 (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 × 1010 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 = T2016-2020 - T2001-2005).

Short-distance movement of individuals might allow buffering the severity of warming impacts on populations, if suitable climatic conditions occur nearby36. 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 (Tvar). 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 (Tbp) was calculated as: Tbp = Tvar/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 (|Tbp| ≥ 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*2 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 R2 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).