Articles | Volume 18, issue 2
Research article
12 Feb 2024
Research article |  | 12 Feb 2024

Recent warming trends of the Greenland ice sheet documented by historical firn and ice temperature observations and machine learning

Baptiste Vandecrux, Robert S. Fausto, Jason E. Box, Federico Covi, Regine Hock, Åsa K. Rennermalm, Achim Heilig, Jakob Abermann, Dirk van As, Elisa Bjerre, Xavier Fettweis, Paul C. J. P. Smeets, Peter Kuipers Munneke, Michiel R. van den Broeke, Max Brils, Peter L. Langen, Ruth Mottram, and Andreas P. Ahlstrøm

Surface melt on the Greenland ice sheet has been increasing in intensity and extent over the last decades due to Arctic atmospheric warming. Surface melt depends on the surface energy balance, which includes the atmospheric forcing but also the thermal budget of the snow, firn and ice near the ice sheet surface. The temperature of the ice sheet subsurface has been used as an indicator of the thermal state of the ice sheet's surface. Here, we present a compilation of 4612 measurements of firn and ice temperature at 10 m below the surface (T10 m) across the ice sheet, spanning from 1912 to 2022. The measurements are either instantaneous or monthly averages. We train an artificial neural network model (ANN) on 4597 of these point observations, weighted by their relative representativity, and use it to reconstruct T10 m over the entire Greenland ice sheet for the period 1950–2022 at a monthly timescale. We use 10-year averages and mean annual values of air temperature and snowfall from the ERA5 reanalysis dataset as model input. The ANN indicates a Greenland-wide positive trend of T10 m at 0.2 C per decade during the 1950–2022 period, with a cooling during 1950–1985 (−0.4C per decade) followed by a warming during 1985–2022 (+0.7 per decade). Regional climate models HIRHAM5, RACMO2.3p2 and MARv3.12 show mixed results compared to the observational T10 m dataset, with mean differences ranging from −0.4C (HIRHAM) to 1.2 C (MAR) and root mean squared differences ranging from 2.8 C (HIRHAM) to 4.7 C (MAR). The observation-based ANN also reveals an underestimation of the subsurface warming trends in climate models for the bare-ice and dry-snow areas. The subsurface warming brings the Greenland ice sheet surface closer to the melting point, reducing the amount of energy input required for melting. Our compilation documents the response of the ice sheet subsurface to atmospheric warming and will enable further improvements of models used for ice sheet mass loss assessment and reduce the uncertainty in projections.

1 Introduction

The Arctic is warming more than 4 times as fast as the global average (Chylek et al., 2022; Rantanen et al., 2022). Consequently, the Greenland ice sheet is exposed to an increase in air temperature (e.g., Hanna et al., 2021; Zhang et al., 2022) which is concurrent to an increase in anticyclonic, cloud-free conditions in summer (Hofer et al., 2017; Ryan et al., 2022). In the low-elevation bare-ice area of the ice sheet, the warming atmosphere increases the sensible heat transfer to the surface (e.g., Wang et al., 2021), while a reduction in cloud cover increases the downward shortwave radiation, both having resulted in melt increases since the late 1980s (Hofer et al., 2017; Trusel et al., 2018; Ryan et al., 2022). Enhanced melt in the bare-ice area initiates multiple feedback processes, such as snowline retreat (Noël et al., 2019; Ryan et al., 2019) and algal growth (e.g., Stibal et al., 2017; Cook et al., 2020), which lead to further expansion and darkening of the bare-ice area and enhanced shortwave radiation absorption. At higher elevations, increased surface melt also triggers a melt–albedo feedback through which liquid water within snow and grain coarsening decreases the snow albedo and increases the absorption of solar radiation (e.g., Nolin and Stroeve, 1997; Box et al., 2012). The increase in ice sheet surface energy influx leads to an increase in surface melt but also to an increase in subsurface temperatures through heat conduction and refreezing of meltwater (Humphrey et al., 2012; Polashenski et al., 2014; McGrath et al., 2013). The subsurface temperature is therefore a key indicator of how the Greenland ice sheet has been affected by recent climatic changes. Furthermore, ice sheet subsurface warming brings the near-surface snow and firn (multi-year, compressed snow) closer to the melting point and makes them less efficient at refreezing and retaining meltwater (Pfeffer et al., 1991; Vandecrux et al., 2020a). Subsurface warming could also trigger thermal-regime shifts across the ice sheet (Marshall, 2021) and increase the ice viscosity (Phillips et al., 2010, 2013; Colgan et al., 2015), although with limited impact on dynamic mass loss (Poinar et al., 2017).

Over the last century, research teams have reported snow, ice and firn subsurface temperatures of the Greenland ice sheet. Of all depths measured, here we focus on measurements at or close to the 10 m depth. The temperature at this depth has been shown to be less affected by seasonal temperature variation and more representative of the long-term temperature and snowfall history at the surface (Dahl-Jensen et al., 1998; McGrath et al., 2013; Kjær et al., 2021). This makes it a convenient standard depth to compare temperatures from different periods and different sites. Here, we compile a dataset of 4612 observations of firm and ice temperature at 10 m below the surface (T10 m), spanning from 1912 to 2022, from published and unpublished sources. We then use 4597 observations of T10 m within the current ice sheet extent and the period 1950–2022 to train an artificial neural network (ANN) model that can predict T10 m over the entire ice sheet. For a given month and location, the ANN estimates T10 m based on 14 parameters derived from the ERA5 reanalysis (Hersbach et al., 2020) that represent the long-term and recent history of air temperature and snowfall. Using our observational dataset of subsurface temperature, as well as our ANN, we evaluate three regional climate models (RCMs) widely used to estimate the surface mass balance of the Greenland ice sheet: RACMO2.3p2 coupled to an offline firn model IMAU-FDM v1.2G (hereafter RACMO, Noël et al., 2019; Brils et al., 2022), MARv3.12 (hereafter MAR, Fettweis et al., 2017, 2020) and HIRHAM5 (hereafter HIRHAM, Langen et al., 2017). We then evaluate the ANN and RCMs' T10 m magnitudes and trends in the bare-ice, percolation and dry-snow areas of the ice sheet. Lastly, we discuss the impact of this subsurface warming on the ice sheet mass balance processes.

2 Methods

2.1 Observed ice sheet subsurface temperature compilation and interpolation

A total of 4612 T10 m observations were compiled from 48 sources (Fig. 1, Table 1). Each dataset is described in the related reference in Table 1, except two yet undescribed datasets. The first unpublished dataset was collected by the late Konrad Steffen and his team and consists of two thermistor strings: one at Swiss Camp, central western Greenland, and another at Summit station, central Greenland, to complement the Greenland Climate Network (GC-Net) automated weather stations (AWSs) at those sites (Steffen et al., 1996; Steffen and Box, 2001). The 11 m long string at Swiss Camp operated between 1992 and 2009 and was equipped with UUB thermistors at 0.5, 0.75 and 1–11 m depth, with 1 m spacing. The 10 to 15 m long string at Summit was equipped with Campbell Scientific T107 thermistors and was active during the periods 2000–2002 and 2007–2009. New sensors were added to the Summit string over the years. The sensors' depth and surface height evolutions could be recovered from field notes, and these data are now presented for the first time. The second unpublished dataset comes from 14 new AWSs installed in 2021 and 2022 by the Geological Survey of Denmark and Greenland (GEUS) as a continuation of the GC-Net sites (Steffen et al., 1996; Steffen and Box, 2001; Vandecrux et al., 2023b). They are equipped with a GeoPrecision TNode thermistor string with sensors installed at 0.5, 1, 1.5, 2, 2.5, 3, 4, 6, 8 and 10 m depth. These data are hosted on the same dataset as the PROMICE AWS data (How et al., 2022).

We also post-processed two previously published datasets. The data from Humphrey et al. (2012) were corrected for the changing depth of the sensor below the surface as snow accumulates or melts at the surface (Sect. S1 in the Supplement) – similarly to the processing of the other time series. The FirnCover dataset (MacFerrin et al., 2022) appeared to have a warm bias due to the use of uncalibrated resistance temperature detectors instead of the conventional thermistor or thermocouple instruments. Using firn temperature observations reported by Samimi et al. (2021) and Heilig et al. (2018) at DYE-2 as a reference, we built an ad hoc correction function that was then applied at all sites within the FirnCover dataset. The correction procedure is described in Sect. S2 and reduces the FirnCover temperatures by 1.1 C on average.

For the temperatures continuously recorded by thermistor or thermocouple strings, the depths of each temperature sensor below the surface were calculated using installation depths and recorded surface height. Wherever necessary, we interpolated the available temperature profiles linearly to 10 m depth and allowed linear extrapolation if at least two measurements were available within 2 m of the 10 m depth. The resulting T10 m values were then aggregated as monthly means if they originated from continuous measurements, or they were left as instantaneous values otherwise.

The measurements conducted by different scientific teams at the same location allow for an assessment of the uncertainty and reproducibility of local vertically interpolated T10 m observations. From 10 sites where simultaneous measurements are available, the median root mean square difference (RMSD) is 0.5 C (Table S1 in the Supplement). Among these 4612 T10 m observations, 15 measurements are either outside of the current ice sheet extent as defined by the GIMP ice sheet delineation (Howat et al., 2014) or outside of the 1950–2022 period we consider for our T10 m reconstruction. There are therefore 4597 T10 m observations in our compilation that can be used for the reconstruction of T10 m on the ice sheet between 1950 and 2022.

Figure 1Spatial (a) and temporal (b) distribution of the T10 m observations in Greenland. Greenland surface classification according to Vandecrux et al. (2019) based on firn density profiles and remote sensing observations.

Table 1Overview of T10 m datasets used in this study.

* Monthly mean values derived from the original measurements

Download Print Version | Download XLSX

2.2 The artificial neural network

Point observations of T10 m only give a partial description of the subsurface temperature: they are discontinuous in space and time. To describe the evolution of T10 m over the entire ice sheet and over the last decades, one can train a machine learning model that links T10 m to an input dataset which is itself continuous in space and time and assumed to drive changes in T10 m. Once the relationship between input and T10 m is learned by the algorithm, the algorithm can be driven by the entire input dataset to reconstruct the T10 m, even at places where no observations are available.

Among machine learning algorithms, ANNs have proven their ability to learn non-linear relationships between a target variable and a set of input variables in numerous glaciological and meteorological applications (e.g., Steiner et al., 2005; Braakmann-Folgmann and Donlon, 2019; Xu et al., 2021). Given that our T10 m compilation, which will be used to train the algorithm, does not encompass all possible ice sheet conditions, we favor ANNs over tree-based algorithms that can be limited when used beyond their training dataset (e.g., Xiong et al., 2020; Liu et al., 2022). Finally, during our search for the most straightforward ANN structure capable of modeling our dataset, we ultimately chose a multi-layer perceptron (Rumelhart et al., 1986). We want to highlight that the choice of model structure, input parameters and training strategy does not have a single optimal configuration. Some of our choices are even decreasing the apparent performance of the model in order to avoid overfitting and to increase the model's capacity to extrapolate outside of its training set. The following sections describe selection of inputs; our enhancement of the dataset's representativity; and eventually the ANN structure, training and uncertainty assessment.

2.2.1 The input parameters

Our target variable, T10 m, is predominantly controlled by (1) the surface temperature through molecular heat conduction, (2) the subsurface refreezing of meltwater through latent heat release and (3) snowfall rates which determine the vertical advection velocity in the firn column. The near-surface air temperature can act as a proxy for both surface temperature and surface melt in the absence of reliable estimates because they all interact within the surface energy budget. The surface temperature itself depends on the near-surface air temperature through turbulent heat fluxes and the surface energy budget. This relationship between air temperature, snowfall and T10 m is notably non-linear. In regions where surface melt is common, meltwater refreezing at depth will lead to T10 m several degrees higher than the average air temperature (e.g., Humphrey et al., 2012). On the other hand, during periods of minimal or no melting (wintertime or nighttime in the summer), the radiative imbalance at the surface and the presence of a near-surface atmospheric temperature inversion can cause the surface temperatures, and, through conduction, the T10 m, to be several degrees lower than the near-surface air temperature (e.g., Miller et al., 2017; Steffen and Box, 2001). Additionally, snowfall affects the subsurface temperature in several ways. In the ablation area, the seasonal snowpack insulates the underlying ice. In the accumulation area, snow accumulated at the surface is, after some time, advected to greater depth, where it can act as either a heat source or sink depending on its temperature at the time of deposition.

Here, we use the air temperature and snowfall monthly grids from the ERA5 reanalysis (Hersbach et al., 2020) to derive our 14 input parameters. We use ERA5 Land at a spatial resolution 0.1×0.1 for 1950–2022 (Muñoz Sabater, 2019) and the original ERA5 (Hersbach et al., 2023) at 0.25×0.25 resolution resampled linearly to 0.1×0.1 for 1940–1950. Delhasse et al. (2020) showed that daily ERA5 near-surface air temperatures compare well with measurements from ice sheet weather stations (mean bias of 0.01 C, root mean square error of 3.05 C). Loeb et al. (2022) found that ERA5's precipitation had the best performance out of three evaluated reanalysis datasets compared to weather station observations in the Canadian Arctic and in Greenland. Using airborne radar measurements of snow accumulation, Ryan et al. (2020) found that ERA5's annual snowfall in Greenland was comparable to estimates from state-of-the-art RCMs and outperformed satellite estimations.

The 10-year average temperature (T2m,10yr) and snowfall (SF10yr) were calculated for each cell and each month to represent the long-term conditions at a given time and place. Additionally, for each grid cell and monthly time step, we calculate the amplitude of the 2 m air temperature in the previous year (T2 m, amp), as well as the average air temperature and snowfall of the 5 previous years. This reflects the capacity of the subsurface to respond not only to long-term changes but also to recent changes in air temperatures and snowfall (e.g., Polashenski et al., 2014). Lastly, to assist the ANN in capturing the annual periodicity, we give as input the cosine of the month (assigning 1 in January and −1 in July). For a given time and location, the ANN therefore takes 14 inputs: T2m,10yr, SF10yr and T2 m, amp; the 5 previous years of annual snowfall; the 5 previous years of air temperature; and the month's cosine.

2.2.2 Weighting of the observations prior to ANN training

Many machine learning algorithms, including ANNs, assume that the training data are representative of the target area (where the model is applied for predictions), i.e., that the data are drawn from the same distribution. This assumption is violated in practice when applying the model to new spatial domains that may contain local conditions not present in the training data. Thus, the representativity of the training dataset compared to the target area is critical for the robustness of any machine learning model, i.e., how well the model generalizes to new and unseen data (Meyer and Pebesma, 2021; Bjerre et al., 2022). The representativity of the 4597 observation sites (training data) compared to the entire Greenland ice sheet where the ANN is applied (target area) was quantified using histogram analysis (Fig. 2). For the three input parameters that define the climate at a given location (T2m,10yr, SF10yr and T2 m, amp), here referred to as pi=1,2,3, we plot the probability histogram of the parameter pi as it appears in ERA5 at our observation locations: this is the observation histogram Ho(pi). We then plot, for that input parameter pi, the probability histogram of all the ice sheet pixels and all time steps within the ERA5 dataset: this is the target histogram Ht(pi). The observation histograms Ho(pi) represent the distribution the ANN will learn from, while the target histograms Ht(pi) represent the values over which the ANN will eventually be applied (Fig. 2). In an ideal scenario where the observational dataset is representative of the parameter space where the ANN will be applied, Ho(pi) and Ht(pi) should show similar distributions.

In practice, the available observations are not representative for the entire ice sheet, stemming from, e.g., monitoring sites producing data continuously or western Greenland being more accessible than eastern Greenland. To make the training dataset more representative of the parameter space in which the ANN will be used, we define for each observation a weight wobs as follows. For each observation and for a given input parameter pi, wobs(pi) is equal to the ratio of Ht(pi) and Ho(pi) for the bin where the observation is located. Consequently, if, in a given bin, the observation histogram is lower than the target histogram, meaning that this bin is underrepresented in the observational dataset compared to the target space, then the weight wobs(pi) will be greater than 1. Inversely, the weight wobs(pi) will be less than 1 if the observation histogram is greater than target histogram. Eventually, for each observation, we calculate the overall weight wobs as the mean of wobs(p1),wobs(p2) and wobs(p3). This overall weight wobs for each observation is used to calculate the loss function (in our case the mean squared error) minimized during the training of the ANN. As a consequence, observations that are located in underrepresented regions of the parameter space will have overall weights wobs>1 and will be given more importance in the training of the ANN. Inversely, observations located in underrepresented parts of the parameter space will have overall weights wobs<1 and will count less in the training of the ANN.

As an illustration, let us consider a T10 m observation from a site and time that has T2m,10yr=-28C. Figure 2a indicates that only ∼10 % of our observation sites have such an average temperature compared to ∼23 % of the ice sheet pixels in ERA5; i.e., this sample comes from an underrepresented temperature range. Following our procedure, we allocate to this observation wobs(p1)=0.23/0.1=2.3 to increase its final weight wobs, which also considers the observation's representativity with regard to SF10yr and T2 m, amp. Inversely, 25 % of our observations have T2m,10yr=-18C, while only 10 % of the ice sheet (according to ERA5) has such an average temperature (Fig. 2a). Consequently, an observation having such T2m,10yr will receive a wobs(p1)=0.1/0.25=0.4 and will weigh less in the training of our ANN.

To verify that our weighting procedure increases the similarity between Ho(pi) and Ht(pi), we evaluate the distance between two histograms H1 and H2 calculated on the same n bins with the Canberra distance (Lance and Williams, 1966; Emran and Ye, 2002):


where Hi=1,2k is the value of histogram Hi=1,2 at bin k. The smaller the Canberra distance dCanberra(Ho(pi),Ht(pi)), the more Ho(pi) and Ht(pi) are similar. The Canberra distance between the observational and target histograms for T2m,10yr,SF10yr and T2 m, amp decreased from 22.7, 12.2 and 14.3 to 11.1, 7.5 and 5.3 when weighing the observations based on their representativity (Fig. 2). Another confirmation that the weights increase the similarity between the observation and target histograms is the smaller difference between the observation and target distributions' median values once the weights are applied: from 4.9 C, 6.0 mm w.e. and 2.4 C with equal weights (Fig. 2a–c) to 2.1 C, 1.4 mm w.e. and 0.4 C with weights (Fig. 2d-f) for T2m,10yr,SF10yr and T2 m, amp respectively.

Figure 2Histograms of the input parameters: 10-year average 2 m air temperature, 10-year average snowfall and annual amplitude of monthly 2 m air temperature. The blue histograms are parameter values as they appear at the observation sites, meaning the training data for the ANN, either with all observations weighing the same (a, b, c) or with weights assigned to each observation based on its representativity (d, e, f). The orange histograms are parameter values as they appear in the ice sheet pixels of ERA5, meaning the target data for the ANN. For each pair of target and observation histograms, we calculate the Canberra distance (dCanberra) as a measure of similarity.


2.2.3 ANN structure and training

Multilayer perceptron ANNs are typically composed of an input layer, with as many nodes as input variables; multiple hidden layers containing several nodes; and an output layer. Each node in the hidden layers (i) makes the weighted sum of the outputs of all nodes from the preceding layer and adds a node-specific bias; (ii) applies a simple, layer-specific activation function to the result; and (iii) passes the output of the activation function to all the nodes of the next layer and so forth. During the training of the model, all the weights and biases from all the nodes are optimized to minimize a loss function. This is done iteratively by (i) passing part of the training set through the ANN, (ii) evaluating the difference between the ANN output and the expected result using the loss function, and (iii) updating the weights and biases to reduce the error in the next iteration (a.k.a. backpropagation). This general ANN structure can be adapted in many ways to the dataset and problem it is applied to. Here, we adjust four of the most important hyperparameters of the ANN: the batch size, i.e., which fraction of the sample is given to the ANN for every training cycle; the number of epochs (or training cycles); the number of layers; and the numbers of nodes within those layers. We use the Adam optimizer (Kingma and Ba, 2014), rectified linear unit activation function (0 if the input is below 0 and f(x)=x if the input x is above 0) and mean squared error as loss functions. Those three settings have been used widely in regression problems (Braakmann-Folgmann and Donlon, 2019; Liu et al., 2022; Lorentzen et al., 2022).

We set the hyperparameters of our ANN in three steps. First, we define a validation set made of 633 observations (14 % of the training dataset) from four sites representing different areas of the ice sheet: NASA-E for the dry-snow area, NASA-SE for the percolation area, and Swiss Camp and KAN_M for the bare-ice area; we use these data as a validation set. Secondly, we train an ensemble of ANNs with two layers of 32 nodes each with batch sizes varying from 100 to 5000 (18 irregular steps) and with between 10 and 1000 epochs (8 irregular steps). Each of the 144 settings are run 10 times to account for the stochastic processes within model training, resulting in a total of 1440 ANNs. We assess the average learning curve for each setting: the mean difference (MD) and root mean squared difference (RMSD) of the trained ANN on the training and validation data as a function of epoch numbers (Fig. S1 in the Supplement). We conclude that (i) small batch sizes (<1000) lead to unstable learning curves (Fig. S1a–d), and (ii) large batch sizes (e.g., 5000) cause slightly slower convergence and similar results as batch sizes of 3000 and 4000 (Fig. S1i–n). From our analysis and as a compromise between stability, rapidity of convergence and potential overfitting, we use a batch size of 4000 over 150 epochs for all ANNs trained henceforth. In the third and last step of our hyperparameter tuning, we use the optimal batch size and number of epochs to train 180 ANNs with one, two and three layers of 8, 16, 32, 64, 128 and 256 nodes each (all layers with same number of nodes, each setting repeated 10 times). We see clear improvements (lower RMSD) when moving from a single layer to two layers and from 8 nodes to 64 nodes (Fig. S2). The improvement moving from two to three layers and from 64 to 128 or 256 nodes is marginal and within the stochastic uncertainty (overlapping standard deviations in Fig. S2c–f). To keep the model design as simple as possible, we henceforth use two layers of 64 nodes each.

Additionally, a Gaussian noise layer that adds random noise to the observations is added after the input layer to further prevent overfitting (e.g., An, 1996). Note that both the addition of Gaussian noise and the assignment of weights to the observations will tend to decrease the apparent performance of the ANN (e.g., MD or RMSD from the non-weighted observational dataset) but will produce a more robust output and prevent overfitting. Considering the limited number of observations relative to the target area, the entire Greenland ice sheet, we train our best model using all the available observations weighted according to their representativity. Consequently, there is no hold out or unseen data for model validation. Alternatively, we use a spatial cross-validation approach to measure the performance and uncertainty of the ANN.

2.2.4 Uncertainty estimation of the ANN with spatial cross-validation

Spatial cross-validation is considered the best-practice approach for evaluating the uncertainty of ANNs when dealing with spatial data (e.g., Brenning, 2012). For this purpose, we separated the Greenland ice sheet into 10 regions (Fig. 3c) after Zwally et al. (2012). Each of the 10 regions contain between 95 and 1298 observations, corresponding to 2 % and 28 % of all observations. For 10 iterations, we hold out the observations located in a different region and train an ANN on the remaining observations. We save these 10 models, and, for any new set of input parameters, we use the standard deviation of the 10 models' predictions as a measure of the uncertainty. This uncertainty is never allowed to be below 0.5 C, which is the measurement uncertainty derived in Sect. 2.1. The monthly grids of ANN uncertainty are provided along with our best estimation of T10 m, which is produced by an ANN trained on all available observations.

For a fair evaluation of our ANN against our observational dataset, we first compare our best ANN model, trained on all T10 m observations, to these same T10 m observations. This evaluation does not show how the model would perform on new, unseen data or regions and consequently leads to an overestimation of the ANN performance. We then compare each T10 m observation to the corresponding T10 m predicted by the one cross-validation model that did not use this observation for training. This second evaluation illustrates how the cross-validation ANNs perform on data that were not included in the training set. It contrasts with the first assessment because it evaluates models that were trained only on part of the observation dataset, and it is therefore a conservative estimate of the performance of the best model trained on all T10 m observations.

2.3 Regional climate models

We evaluate 10 m subsurface temperatures from three regional climate models: MARv3.12 (Fettweis et al., 2017, 2020), RACMO2.3p2 (Noël et al., 2019) with the updated IMAU-FDMv1.2G (Brils et al., 2022) and HIRHAM5 (Langen et al., 2017). We calculate the MD and RMSD between the observed and simulated 10 m subsurface temperatures. For this study, the outputs from MAR, RACMO and HIRHAM are available over the periods 1950–2020, 1958–2020 and 1980–2016, respectively. We compare each model to the measurements within the common 1980–2016 period for which all three model outputs are available, as well as against all observations.

All three models use a multilayer snow, firn and ice model to calculate subsurface temperature. In addition to differences in surface forcing in the three models (e.g., in snowfall, rainfall, melt and energy fluxes), the models also differ in the way they calculate the subsurface characteristics that impact the subsurface temperature. Both MAR and HIRHAM estimate firn densification using the overburden pressure, respectively from Brun et al. (1989) and Vionnet et al. (2012), while RACMO uses a compaction law that was derived for steady-state firn (Arthern et al., 2010) and empirically fitted to observations (Ligtenberg et al., 2011; Brils et al., 2022). RACMO's offline run with IMAU-FDMv1.2G uses the thermal conductivity parameterization from Calonne et al. (2019), while HIRHAM and MAR use the parameterization by Yen (1981). The three models treat the release of latent heat during the refreezing of meltwater in a similar manner, but the meltwater infiltration is calculated differently. Both MAR and RACMO use a bucket scheme: meltwater infiltrates downward unless the water is refrozen or retained through capillary forces, and ice layers are considered permeable at the model scale (Ligtenberg et al., 2018). In HIRHAM, the use of a parameterization of Darcy flow (Hirashima et al., 2010) and accounting for the decrease of the layer permeability due to ice content (Colbeck, 1975) lead to shallower infiltration than in RACMO (Vandecrux et al., 2020b). Another model detail that impacts the calculated subsurface temperature is the boundary condition at the bottom of the model domain. HIRHAM uses a temperature scheme that requires a fixed temperature at the lowermost firn layer which is set, for each pixel, to the long-term air temperature average (Langen et al., 2017). Both MAR and RACMO use the Neumann boundary condition at the bottom layer of the firn model, which implies no heat flux through the lower boundary of the model. In the ablation area, new material needs to be provided to the bottom layer of the model as surface ablation melts ice away. In MAR, as soon as the model column height is lower than 29 m, a 1 m thick layer composed of ice is added at the bottom of the model column. MAR then uses a simple assumption that the underlying ice would always be cooler than the ablating, near-surface ice. Consequently, the temperature of the 1 m layer added at the bottom of the model was fixed to be 1 % lower (when calculated in Kelvin) than the temperature of the lowermost layer left in the model. The differences between RCM-simulated subsurface temperatures are partly due to these different modeling approaches for the subsurface processes. This can be illustrated when different subsurface models are forced with similar surface data (Lundin et al., 2017; Vandecrux et al., 2020b). Another source of discrepancy is the difference in surface climate that is simulated in each of these three models. More information about the accuracy of the simulated surface climate within each RCM can be found in Fettweis et al. (2020), Langen et al. (2017) and Noël et al. (2019).

Figure 3(a) Evaluation of theT10 m simulated by the best ANN model against the observations used for training. (b) Evaluation of the T10 m simulated by the 10 cross-validation ANN models against their unseen data (i.e., not used for training). The statistics presented are mean difference (MD), root mean square difference (RMSD) and number of samples for which the comparison was possible (N). (c) The 1950–2022 average of the ANN uncertainty as calculated from the standard deviation of 10 cross-validation ANN models trained on different spatial subsets of the observation dataset. (d–h) Examples of ANN T10 m prediction, its uncertainty and the prediction of the 10 cross-validation models at a dry-snow site (NASA-E), one percolation site (DYE-2), two ablation sites (KAN_L, KPC_U) and a firn aquifer site (FA_13).

Table 2Trends in 10 m subsurface temperature (T10 m) calculated from the ANN and observations (obs.) at 10 sites for the periods 1998–2010 and 1998–2022. ANN trends are calculated only from the months where observations are also available. The difference between the two calculated trends and the number of monthly values used for the calculation (N) are also given for each site.

Download Print Version | Download XLSX

3 Results

3.1 Performance of the ANN

When comparing the best ANN model to the T10 m observations it was trained on, we find an MD of 0.0 C and an RMSD of 1.6 C (Fig. 3a). However, when evaluating the cross-validation models against their respective unseen data, we find a similar MD (0.1 C) and an RMSD of 2.5 C (Fig. 3b). While the first evaluation is overoptimistic, the second does not directly evaluate our best ANN model, which is trained on all available data. These estimates nevertheless provide bounds to the true performance of our ANN.

Averaging over the entire period 1950–2022, the ANN uncertainty is lowest across the dry snow area (Fig. 3c), illustrated by the NASA-E site (Fig. 3d). The uncertainty increases towards the ice sheet margin in west, north and northeast Greenland (Fig. 3c), as exemplified by the site DYE-2 in the percolation area in western Greenland (Fig. 3e) and by the sites KAN_L and KPC_U (Fig. 3f, g), two ablation-area sites in western and northeastern Greenland, respectively. The ANN uncertainty peaks in southeast Greenland (Fig. 3c) where relatively high air temperatures and snow accumulation produce temperate firn conditions and firn aquifers (Forster et al., 2013; Kuipers Munneke et al., 2014). When the measurements in this region are removed for cross-validation, there are no firn aquifer observations left in the training set for the ANN to learn what the T10 m structure in this ice sheet region is. This is illustrated in Fig. 3h by a cross-validation model predicting lower T10 m than observed at the site FA_13, resulting in a larger standard deviation between the cross-validation models for FA_13. Our uncertainty estimation is conservative because the final ANN model is eventually trained on all observations. The regions of high uncertainty highlight where observations are the most needed to map the ice sheet subsurface temperature.

Figure 4Evaluation of the monthly 10 m subsurface temperatures simulated by RACMO (a), HIRHAM (b) and MAR (c) against observations. The statistics presented are mean deviation (MD), root mean square difference (RMSD) and number of samples for which the comparison was possible (N) for the period when all three models are available (1980–2016). For RACMO and MAR, the statistics for all available measurements are given in the parentheses. For sites where annual surface ablation is larger than snow accumulation, i.e., net ablation sites with a bare-ice cover in summer, symbols and statistics are shown in red.


To evaluate the capacity of the ANN to capture the recent evolution of T10 m, we select 10 sites where more than 60 monthly values are available between 1998 and 2022 and compare the trends calculated from the ANN and the observations over the periods 1998–2010 and 1998–2022 (Table 2). These periods were chosen because of a general lack of measurements between 2011 and 2020. Trends calculated from the ANN only consider the months where observations are available. We note that, due to the missing months, these trends are not reliable for general inference of the true T10 m evolution: depending on which months are missing, it might overestimate or underestimate the true T10 m trend for these periods. The median T10 m trends for 1998–2010 are 0.9 and 0.8 C per decade for the ANN and for the observations, respectively (Table 2). For the period 1998–2022, the median T10 m trends for 1998–2010 are 0.4 and 0.6 C per decade for the ANN and for the observations, respectively (Table 2). The ANN therefore slightly overestimates the T10 m trend during 1998–2010 and underestimates it during 1998–2022. We conclude that the ANN reproduces the magnitude of the T10 m increase seen in observations, although this aptitude varies with the location and the time period considered. From this assessment and because the ANN does not suffer either temporal or spatial gaps, the ANN appears to be a suitable tool to study the trends in T10 m over the entire Greenland ice sheet.

3.2 RCM evaluation and comparison with the ANN

We evaluate the RCMs against the observed T10 m in the period 1980–2016, for which all three RCMs' outputs are available (Fig. 4). HIRHAM shows the best performance (MD =−0.4C, RMSD = 2.8 C), followed by RACMO (MD =−1.3C, RMSD = 3.1 C) and MAR (MD =+1.2C, RMSD = 4.7 C). For the observation sites located in the ablation area, RACMO, HIRHAM and MAR have a cold bias with MD of −3.6, −0.9 and −3.4C, respectively (Fig. 4). MAR captures neither the geographical nor the seasonal variability of T10 m in the ablation area (RMSD = 5.4 C). The ANN, although of a different nature, gives better statistics at these ablation sites with an MD of 0.2 C and an RMSD of 2.9 C, even when calculated from our cross-validation models' unseen data (Fig. 3b).

We further evaluate the ANN and RCMs at eight sites (Fig. 5) that are representative of the dry-snow (Summit, NASA-E), percolation (DYE-2, KAN_U), bare-ice (Swiss Camp, KPC_U, SCO_U) and firn aquifer regions (FA_13). The ANN performs well at most of these sites: the average MD for these eight sites is less than 0.2 C, and the average RMSD is 1.2 C. RACMO overestimates T10 m at lower-temperature sites in the dry-snow area (Fig. 5a, b) and underestimates T10 m at the accumulation sites with relatively high melt (Fig. 5c, d) and at ablation sites (Fig. 5e, g). HIRHAM compares better than RACMO to the measurements at accumulation sites (Fig. 5a, b) and can either over- or underestimate T10 m at percolation sites (Fig. 5c, d) and at ablation sites (Fig. 5e–g). MAR simulates T10 m that are unrealistic both in magnitude and in variation (Fig. 5). The causes of this low performance will be discussed in Sect. 4. At a firn aquifer site (Fig. 5h), the ANN and the three RCMs successfully estimate relatively high T10 m during the period 2013–2015 for which observations are available. Yet, the models diverge significantly when estimating the past history of the site: HIRHAM and MAR indicate T10 m close to 0 C from the models' respective initiations in 1980 and 1950, while RACMO and the ANN indicate that T10 m below −2C may have been common at FA_13 before 2000 (Fig. 5h).

Figure 5Observed and simulated 10 m subsurface temperatures at selected sites. Note the different y axes.


3.3T10 m trends in the ANN and RCMs

According to the ANN, the Greenland ice sheet average T10 m has been increasing significantly at a rate of +0.2C per decade (P<0.01) over the 1950–2022 period (Table 3, Fig. 6a), from an ice-sheet-wide average value of −21.1C in 1950 to −19.2C in 2022. This increase was not constant over the 1950–2022 period. When fitting multiple piecewise linear functions to the Greenland ice sheet average T10 m, with a breakpoint between 1951 and 2021, we identify 1985 as the breakpoint year that explains most of the variance in the ice-sheet-wide average T10 m time series. This piecewise linear function consists of a period of significant cooling between 1950 and 1985 (−0.4C per decade, P<0.01) followed by a strong warming from 1985 to 2022 (+0.7C per decade, P<0.01). Both the cooling that occurred until 1985 and the subsequent warming were most pronounced in central and southern Greenland (Fig. 7a, b). In contrast, the low elevations of the northwest Greenland ice sheet underwent warming during the entire period (Fig. 7a–c).

For the time period for which ANN, RACMO, HIRHAM and MAR are available (1980–2016), the ANN gives an ice-sheet-wide average T10 m trend of +0.6C per decade (P<0.01), while the equivalent trends are estimated at +0.3, +0.4 and −0.1C per decade by RACMO, HIRHAM and MAR (P≤0.01), respectively (Table 3, Fig. 6a). The spatial patterns of T10 m trends in the three RCMs (Fig. 7e–g) are consistent with the ANN (Fig. 7d): a more pronounced warming at a mid-elevation band around the ice sheet and a milder warming (or cooling for MAR) in the rest of the ice sheet.

Since the processes controlling T10 m depend on the local climatic, snow and ice conditions, we also compare the evolution of T10 m in different ice sheet regions (Fig. 1): (i) the bare-ice area where seasonal snow melts completely and exposes underlying glacial ice at the end of summer, (ii) the dry-snow area where little or no melt occurs, and (iii) the intermediate percolation area where a significant portion of the annual snow accumulation melts in spring and summer and percolates into the underlying firn (Fig. 1a). In the bare-ice area (Fig. 6b), the observation-based ANN predicts stable T10 m until the 1980s and increasing T10 m thereafter (+0.6C per decade over 1985–2022, P<0.01). In contrast, MAR estimates a negative trend in T10 m temperatures over the 1950–2022 period and overestimates the T10 m during the 1950–2000 period compared to the ANN (Fig. 6b). In the bare-ice area, RACMO and HIRHAM both present a T10 m trend of +0.2C per decade (P<0.01) over 1980–2016 period, which is 66 % smaller than the ANN trend in the ablation area for the same period. In the dry-snow area (Fig. 6c), there is a better agreement between the models, but lower T10 m in the 1990s in the ANN leads to a more pronounced warming trend in that area (+0.5C per decade over 1980–2016, P<0.01), which is 40 %–60 % larger than warming trends predicted by RACMO and HIRHAM. MAR describes an overall cooling in the dry-snow area (Figs. 6c, 7g, Table 3). In the percolation area (Fig. 6d), MAR has a warm bias compared to the other models (+4C on 1980–2016), but all models agree on the strong warming that has occurred here since the 1980s: between +0.5 and +0.9C per decade (all P<0.01) over 1980–2016 (Figs. 6d, 7d–g, Table 3). Overall, these spatial differences average into a warm bias of MAR for the entire ice sheet and more pronounced trends for the ANN than for the RCMs (Fig. 6a).

Figure 6Evolution of the 10 m subsurface temperature (T10 m) for all of the Greenland ice sheet (a) and in three ice sheet regions (b–d). Although all panels have the same vertical-axis scaling, note the different vertical-axis bounds.


Figure 7Trends in 10 m subsurface temperature as determined by the ANN over the periods 1950–1984 (a), 1985–2022 (b), 1950–2022 (c) and 1980–2016 (d) and calculated by RACMO (e), HIRHAM (f) and MAR (g) over the period 1980–2016, when data from all models are available. Dotted areas indicate trends below significance level (P>0.1). In panel (b), the lower limits of the dry-snow area (DSA) and of the percolation area (PA) are shown in dark green and black, respectively.

Table 3Trends in 10 m subsurface temperature for different ice sheet regions and different periods. All trends are significant at a 0.1 level.

Download Print Version | Download XLSX

4 Discussion

We compiled the largest dataset of observed subsurface temperature on the Greenland ice sheet to date and used it to train an ANN which, with snowfall and temperature from ERA5 reanalysis as input, estimates monthly grids of 10 m subsurface temperature over the entire ice sheet for the 1950–2022 period. The ANN describes a −0.4C per decade T10 m trend during 1950–1985 (Fig. 6a, Table 3), which is consistent with the negative trends in air temperatures found by Zhang et al. (2022) in the ERA5 reanalysis and RACMO RCM from the late 1950s to early 1990s. The following increase in T10 m (+0.7C per decade) calculated by the ANN from the 1990s to 2022 is consistent with all-year and summer air temperature increases found in weather station observations, reanalysis datasets and regional climate models (Hanna et al., 2021; Zhang et al., 2022). Between 1950 and 2022, the average T10m trend derived from our observation-based ANN is identical to the trend in annual air temperature in ERA5, showing that the subsurface likely responds to the increasing air temperature at the surface. Additionally, the ANN estimates a strong warming of +0.9C per decade on average and up to +1.4C per decade locally in the percolation area (bounded by the dark-green and black lines in Fig. 7b) during the 1985–2022 period.. This localized warming of the percolation area is also calculated by the three RCMs (Figs. 6d, 7e–g). However, this hotspot of T10 m increase is not found in air temperature trends (Zhang et al., 2022, Figs. 5–7). The warming of the subsurface in the percolation area stems from the increased meltwater infiltration and from the latent heat released by refreezing (e.g., Humphrey et al., 2012; Vandecrux et al., 2020a). This successful identification of areas subject to firn warming by the ANN is remarkable considering that the ANN only learns from the T10 m observations and the local air temperature and snowfall history and is not fed information on meltwater infiltration and refreezing. This indicates that the ANN successfully learns which areas are susceptible to undergoing meltwater infiltration and refreezing from its training data.

The ANN model has the strength of statistical models: it fits the training data and thereby performs better than RCMs when evaluated against observations used for training (Figs. 3a, 4). Yet, ANNs and statistical models have several limitations. Firstly, the ANN is greatly dependent on the distribution of the training data and how representative those data are of the parameter space where the ANN is applied. Our methods give more weight to observations that are from underrepresented areas of the parameter space. Yet, there are still regions with particular combinations of air temperature and snowfall where no observations are available and where the ANN extrapolates. More observations are needed from these less-visited parts of the ice sheet to further train the ANN. These new measurements could either focus on the coldest parts of the ice sheet, where our compilation currently lacks measurements (Fig. 2a), or on the areas where our uncertainty is the highest, in the southeast (Fig. 3c). Secondly, the ANN is limited by the input parameters it draws on. For instance, inaccuracies in ERA5 data (as discussed in Delhasse et al., 2020; Zhang et al., 2022) for certain periods or locations will affect the performance of the ANN, as will T10 m measurement uncertainties. Besides, only using two input parameters (air temperature and snowfall) must introduce inaccuracy through oversimplification of complex physical processes. Additionally, the relatively coarse resolution of the input grid (0.1×0.1) prevents the ANN from identifying local phenomena such as localized meltwater refreezing in surface deepenings and crevasses (Hills et al., 2018; Chudley et al., 2021) or in ephemeral perched aquifers (Humphrey et al., 2021; Culberg et al., 2022). Furthermore, our ANN model cannot account for the exposure of ice affected by past temperature anomalies, i.e., the advection of deep ice in the ablation zone that may drive T10 m more than surface conditions (Lüthi et al., 2015). Other widespread processes such as the penetration of shortwave radiation into the subsurface (Van den Broeke et al., 2008; Kuipers Munneke et al., 2009; Van Dalum et al., 2021), firn ventilation (Albert and Shultz, 2002) or potentially wind pumping (Clarke et al., 1987) are more likely to be accounted for by the ANN because observations subject to these processes are included in the dataset. Ultimately, the ANN cannot identify the processes that are responsible for a given subsurface temperature, but it can learn which T10 m are usually seen at various temperature and snowfall combinations.

Although the RCMs calculate subsurface temperatures in similar ways (see Sect. 2.3), differences arise due to their various assumptions. For example, MAR assumes that, in the ablation area, the material added at the lower bound of the model column is always slightly colder than the lowermost material left in the column. This explains the decreasing trend in simulated T10 m in Fig. 6a and the inability of MAR to explain the observed T10 m variation at the ablation sites (Fig. 4c). The noise within the T10 m trend map (Fig. 7g) is also indicative of some numerical instability in this deep temperature prescription. These limitations of the model's boundary conditions have now been identified, and efforts are ongoing to remediate them in the next version of MAR. It is, however, interesting to note that these biases do not significantly impact the surface mass balance simulated by MAR; different sensitivity tests were performed with the aim of improving the comparison with T10 m, and for all of them, the MAR results at the surface remained unchanged. Langen et al. (2017) showed that the simulated subsurface temperature profile in HIRHAM in the percolation area greatly depends on the formation, in the model, of ice layers of density greater than 830 kg m−3 that inhibit water infiltration. The formation of these high-density layers in the model depends on the surface climate and subsurface model but also on the discretization of the modeled firn column, which is currently fixed in HIRHAM (Vandecrux et al., 2020b). Recent efforts to update the HIRHAM subsurface scheme to a more flexible discretization that would preserve high-density layers were made for Antarctica (Hansen et al., 2021) but have not yet been applied to Greenland. Steger et al. (2017) found that the SNOWPACK model forced by RACMO2.3, an older version of RACMO2.3p2, evaluated here, overestimated the subsurface temperature in the high-elevation areas in northwestern Greenland while underestimating the firn temperature at lower elevations due to either insufficient meltwater generation at the surface or too-shallow simulated meltwater infiltration. In that same study, RACMO2.3, in combination with both IMAU-FDMv1.1 and SNOWPACK subsurface schemes, could not accurately reproduce subsurface temperatures at some low-elevation firn sites because the models represented them as bare-ice sites. This mismatch between the simulated and actual surface type – bare ice or porous firn – makes sites at the transition between the bare-ice and percolation areas, i.e., the equilibrium line, particularly challenging for all RCMs (e.g., KAN_U, Swiss Camp, KPC_U in Fig. 5d–f). Switching from version 2.3 to 2.3p2, in combination with an update to IMAU-FDMv1.2, allowed RACMO to simulate KAN_U as a firn site rather than a bare-ice site (Ligtenberg et al., 2018; Brils et al., 2022). The IMAU-FDM always allows meltwater infiltration, which may lead to an overestimation of T10 m at sites where thick ice layers in the firn provide a barrier for further percolation. This was highlighted at KAN_U when driving IMAU-FDMv1.1 with surface temperature and melt rates derived from observations (Vandecrux et al., 2020b). However, the updated IMAU-FDMv1.2 forced by RACMO2.3p2 now shows a slight cold bias at KAN_U (Fig. 5d), indicating that too-deep meltwater infiltration is no longer an issue at that site (Brils et al., 2022).

The subsurface temperature impacts the surface energy budget through the conductive heat flux and thereby affects the snow and ice surface melt. Heat from a warm subsurface will be conducted to the surface when surface temperatures are lower. And vice versa, a colder subsurface represents a heat sink (heat will conduct down, away from the surface) and will moderate surface melt. Another consequence of the near-surface snow and firn warming is that it decreases the cold content and therefore the retention capacity of the snow and firn (Pfeffer et al., 1991; Vandecrux et al., 2020a). Meltwater retention in firn occurs when (i) pore space is available, (ii) this pore space can be accessed by the meltwater, and (iii) cold content is available to refreeze the meltwater. Vandecrux et al. (2019) documented that the upper 10 m of the firn in the lower part of the accumulation area has lost about 20 % of its pore space over the last decades, while, in the dry-snow area, pore space has remained stable since the 1950s. Our work documents the recent subsurface warming of the ice sheet and how the upper 10 m of snow, ice and firn is brought closer to the melting point, potentially enhancing meltwater runoff in the subsequent summers. Our ANN estimates that the dry-snow area average T10 m increased from −27.3C over 1980–1990 to −25.8C over 2010–2020. Similarly, the percolation area warmed from −13.7C over 1980–1990 to −11.8C over 2010–2020, 1.9 C (14 %) closer to the melting point. Our findings complement other work showing the changes in the Greenland ice sheet subsurface and their impact on the ice sheet mass loss: for example, the recent expansion of the firn aquifer area stemming from, among other causes, the loss of firn cold content (Horlings et al., 2022) or the increasing the runoff from the firn area (Tedstone and Machguth, 2022) linked to the formation of ice layers reducing meltwater percolation and retention into the underlying firn (Machguth et al., 2016; MacFerrin et al., 2019) .

5 Conclusion

Using the most complete compilation of observed T10 m on the Greenland ice sheet to date, we trained an artificial neural network (ANN) to describe the spatio-temporal evolution of T10 m during 1950–2022. We found that, following a significant cooling between 1950 and 1985 (−0.4C per decade, P<0.01), ice-sheet-wide T10 m increased by +0.7C per decade from 1985 to 2022 (P<0.01). Overall, the Greenland ice sheet T10 m increased at a rate of +0.2C per decade over the 1950–2022 period in response to increasing energy influx at the surface. Our observational dataset yielded unique and extensive constraints on the subsurface temperature simulated by three conventional regional climate models, RACMO, MAR and HIRHAM, and demonstrated their mixed performance. Notably, it revealed numerical instabilities in MAR, prompting improvements in its snow module, although these T10 m biases apparently have low impact on the SMB simulated by MAR. This work highlights the value of in situ measurements of ice, snow and firn temperatures to better quantify the response of the Greenland ice sheet to Arctic warming and to reduce uncertainty in projections of mass loss from the Greenland ice sheet. Our evaluation shows highest ANN uncertainty in the southeast and in the lower percolation area in northern Greenland (Fig. 3). Those are regions where few observations are available (Fig. 1)and where additional measurements are needed to better understand the ice sheet's subsurface temperature and its response to recent climate change.

Code and data availability

The original subsurface temperature datasets are cited in Table 1, and, when available, download links to the original datasets can be found in the reference list. Some of the observational data were found in other compilations such as Mock and Weeks (1965), McGrath et al. (2013) or Løkkegaard et al. (2023), and Mankoff et al. (2022).

The monthly T10 m dataset is currently hosted at (Vandecrux, 2023a) and has also been added to the 2023 release of the SUMup dataset (, Vandecrux et al., 2023a). A compilation of non-interpolated, instantaneous subsurface temperatures can be requested from the authors. The T10 m maps here are available at (Vandecrux, 2023b), and the scripts used for the analysis are available at (Vandecrux, 2023c).


The supplement related to this article is available online at:

Author contributions

BV designed the study, processed and compiled the observations, and conducted the analysis with help from RSF, JEB, FC, RH and ÅKR. RSF, APA, JEB, FC, RH, ÅKR, AH, JA, PCJPS and DvA funded or organized projects measuring subsurface temperatures and provided valuable data. EB provided insights into the ANN setup and tuning. XF, PKM, MRvdB, MB, PLL and RM provided RCM outputs and their insights on their performance. BV drafted the paper, to which all the co-authors contributed.

Competing interests

At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We salute the field parties that collected data in the field over more than 100 years, allowing us to conduct the present analysis. We thank Hubertus Fischer, Wilfried Häberli, Olaf Eisen, Shin Sugiyama, Sumito Matoba and Robert Hawley for sharing their subsurface temperature data. We also thank Valentina Radić for the editing of our paper and two anonymous reviewers for their constructive comments.

Financial support

This research has been supported by the PROMICE and GC-Net programmes (, last access: 8 February 2024) which are funded by the Danish Ministry for Climate, Energy and Utilities. The historical GC-Net AWS and FirnCover measurements were supported by NASA and NSF grants. Jakob Abermann was supported by the Austrian Science Fund (grant no. P35388). Federico Covi, Åsa K. Rennermalm and Regine Hock were supported by US National Science Foundations (NSF) grant nos. 397516-66782, OPP-1604058 and OPP-1603815, and by The United States Ice Drilling Program through the NSF Cooperative Agreement (grant no. 1836328). Michiel R. van den Broeke, Peter Kuipers Munneke and Max Brils received funding from the Netherlands Earth System Science Centre (NESSC). Peter L. Langen received funding from the Aarhus University Interdisciplinary Centre for Climate Change (iClimate, Aarhus University). Achim Heilig has been supported by the DFG (Deutsche Forschungsgemeinschaft; grant no. HE7501/1-1).

Review statement

This paper was edited by Valentina Radić and reviewed by two anonymous referees.


Abermann, J., Vandecrux, B., Scher, S., Löffler, K., Schalamon, F., Trügler, A., Fausto, R. S., and Schöner, W.: Learning from Alfred Wegener's pioneering field observations in West Greenland after a century of climate change, Sci. Rep., 13, 7583,, 2023. 

Albert, M. R. and Shultz, E. F.: Snow and firn properties and air–snow transport processes at Summit, Greenland, Atmos. Environ., 36, 2789–2797,, 2002. 

Ambach, W.: Zum Wärmehaushalt des Grönländischen Inlandeises: Vergleichende Studie im Akkumulations- und Ablationsgebiet, Polarforschung, 49, 44–54, 1979. 

An, G.: The Effects of Adding Noise During Backpropagation Training on a Generalization Performance, in: Neural Computation, 8, 643–674,, April 1996. 

Arthern, R. J., Vaughan, D. G., Rankin, A. M., Mulvaney, R., and Thomas, E. R.: In situ measurements of Antarctic snow compaction compared with predictions of models, J. Geophys. Res.-Earth Surf., 115, 1–12,, 2010. 

Benson, C. S.: Stratigraphic Studies in the Snow and Firn of the Greenland ice sheet. Research Report 70. Reprinted August 1996 with revision of Cold Regions Research Laboratory (CRRL), 93, (last access: 5 February 2024), 1962. 

Bjerre, E., Fienen, M. N., Schneider, R., Koch, J., and Højberg, A. L.: Assessing spatial transferability of a random forest metamodel for predicting drainage fraction, J. Hydrol., 612, 128177,, 2022. 

Box, J. E., Fettweis, X., Stroeve, J. C., Tedesco, M., Hall, D. K., and Steffen, K.: Greenland ice sheet albedo feedback: thermodynamics and atmospheric drivers, The Cryosphere, 6, 821–839,, 2012. 

Braakmann-Folgmann, A. and Donlon, C.: Estimating snow depth on Arctic sea ice using satellite microwave radiometry and a neural network, The Cryosphere, 13, 2421–2438,, 2019. 

Braithwaite, R.: Firn temperature and meltwater refreezing in the lower accumulation area of the Greenland ice sheet, Pâkitsoq, West Greenland, Rapport Grønlands Geologiske Undersøgelse, 159, 109–114,, 1993. 

Brenning, A.: Spatial cross-validation and bootstrap for the assessment of prediction rules in remote sensing: The R package sperrorest, 2012 IEEE International Geoscience and Remote Sensing Symposium, Munich, Germany, 5372–5375,, 2012. 

Brils, M., Kuipers Munneke, P., van de Berg, W. J., and van den Broeke, M.: Improved representation of the contemporary Greenland ice sheet firn layer by IMAU-FDM v1.2G, Geosci. Model Dev., 15, 7121–7138,, 2022. 

Brun, E., Martin, E., Simon, V., Gendre, C., and Coleou, C.: An energy and mass model of snow cover suitable for operational avalanche forecasting, J. Glaciol., 35, 333–342,, 1989. 

Calonne, N., Milliancourt, L., Burr, A., Philip, A., Martin, C. L., Flin, F., and Geindreau, C.:. Thermal conductivity of snow, firn, and porous ice from 3-D image-based computations, Geophys. Res. Lett., 46, 13079–13089,, 2019. 

Charalampidis, C., Van As, D., Colgan, W. T., Fausto, R. S., Macferrin, M., and Machguth, H.: Thermal tracing of retained meltwater in the lower accumulation area of the Southwestern Greenland ice sheet, Ann. Glaciol., 57, 1–10,, 2016. 

Charalampidis, C., van As, D., Colgan, W., Fausto, R., MacFerrin, M., Machguth, H., and Vandecrux, B.: Subsurface monitoring at 1840 m a.s.l. on the southwestern Greenland ice sheet during 2012–2013, GEUS Dataverse [data set],, 2022. 

Chudley, T. R., Christoffersen, P., Doyle, S. H., Dowling, T. P. F., Law, R., Schoonman, C. M., Bougamont, M., and Hubbard, B.: Controls on water storage and drainage in crevasses on the Greenland Ice Sheet, J. Geophys. Res.-Earth Surf., 126, e2021JF006287,, 2021. 

Chylek, P., Folland, C., Klett, J. D., Wang, M., Hengartner, N., Lesins, G., and Dubey, M. K.: Annual mean Arctic Amplification 1970–2020: Observed and simulated by CMIP6 climate models, Geophys. Res. Lett., 49, e2022GL099371,, 2022. 

Clarke, G. K. C., Fisher, D. A., and Waddington, E. D.: Wind pumping: A potentially significant heat source in ice sheets, Phys. Basis Ice Sheet Model. IAHS Publ., 170, 169–180, 1987. 

Clausen, H. B. and Hammer, C. U.: The Laki and Tambora Eruptions as Revealed in Greenland Ice Cores from 11 Locations, Ann. Glaciol., 10, 16–22,, 1988. 

Clausen, H. B., Gundestrup, N. S., Johnsen, S. J., Bindschadler, R. A., and Zwally, J.: Glaciological investigations in the Crête area, central Greenland: a search for a new deep-drilling site, Ann. Glaciol., 10, 10–15, 1988. 

Clement, P.: Glaciological Activities in the Johan Dahl Land Area, South Greenland, As a Basis for Mapping Hydropower Potential, Rapport Grønlands Geologiske Undersøgelse, 120, 113–21,, 1984. 

Colbeck, S. C.: A theory for water flow through a layered snowpack, Water Resour. Res., 11, 261–266,, 1975. 

Colbeck, S. C. and Gow, A. J.: The Margin of the Greenland ice sheet at Isua, J. Glaciol., 24, 155–165,, 1979. 

Colgan, W. and Vandecrux, B.: Camp Century: Firn temperature measurements (CEN-THM), GEUS Dataverse [data set],, 2021. 

Colgan, W., Sommers, A., Rajaram, H., Abdalati, W., and Frahm, J.: Considering thermal-viscous collapse of the Greenland ice sheet, Earth's Future, 3, 252–267,, 2015. 

Cook, J. M., Tedstone, A. J., Williamson, C., McCutcheon, J., Hodson, A. J., Dayal, A., Skiles, M., Hofer, S., Bryant, R., McAree, O., McGonigle, A., Ryan, J., Anesio, A. M., Irvine-Fynn, T. D. L., Hubbard, A., Hanna, E., Flanner, M., Mayanna, S., Benning, L. G., van As, D., Yallop, M., McQuaid, J. B., Gribbin, T., and Tranter, M.: Glacier algae accelerate melt rates on the south-western Greenland Ice Sheet, The Cryosphere, 14, 309–330,, 2020. 

Covi, F., Hock, R., Rennermalm, A., Leidman, S., Miege, C., Kingslake, J., Xiao, J., MacFerrin, M., and Tedesco, M.: Meteorological and firn temperature data from three weather stations in the percolation zone of southwest Greenland, 2017–2019, Arctic Data Center [data set],, 2022. 

Covi, F., Hock, R., and Reijmer, C.: Challenges in modeling the energy balance and melt in the percolation zone of the Greenland ice sheet, J. Glaciol., 69, 164–178,, 2023. 

Culberg, R., Chu, W., and Schroeder, D. M.: Shallow fracture buffers high elevation runoff in Northwest Greenland, Geophys. Res. Lett., 49, e2022GL101151,, 2022. 

Dahl-Jensen, D., Mosegaard, K., Gundestrup, N., Clow, G. D., Johnsen, S. J., Hansen, A. W., and Balling, N.: Past temperatures directly from the Greenland ice sheet, Science, 282, 268–271,, 1998. 

Davis, R. M.: Approach roads Greenland 1960–1964 Technical Report 133, Corps of Engineers Cold Regions Research and Engineering Laboratory, (last access: 5 February 2024), 1967. 

Davis, T. C.: Structures in the upper snow layers of the southern Dome Greenland ice sheet, CRREL research report 115, (last access: 5 February 2024), 1954. 

Delhasse, A., Kittel, C., Amory, C., Hofer, S., van As, D., S. Fausto, R., and Fettweis, X.: Brief communication: Evaluation of the near-surface climate in ERA5 over the Greenland Ice Sheet, The Cryosphere, 14, 957–965,, 2020. 

de Quervain, M.: Schneckundliche Arbeiten der Internat. Glaziolog. Gronlandexpedition, Meddelelser om Gronl., 177, 283 pp., 1969. 

Emran, S. M. and Ye, N.: Robustness of Chi-square and Canberra distance metrics for computer intrusion detection, Qual. Reliab. Engng. Int., 18, 19–28,, 2002. 

Fausto, R. S., van As, D., Mankoff, K. D., Vandecrux, B., Citterio, M., Ahlstrøm, A. P., Andersen, S. B., Colgan, W., Karlsson, N. B., Kjeldsen, K. K., Korsgaard, N. J., Larsen, S. H., Nielsen, S., Pedersen, A. Ø., Shields, C. L., Solgaard, A. M., and Box, J. E.: Programme for Monitoring of the Greenland Ice Sheet (PROMICE) automatic weather station data, Earth Syst. Sci. Data, 13, 3819–3845,, 2021. 

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033,, 2017. 

Fettweis, X., Hofer, S., Krebs-Kanzow, U., Amory, C., Aoki, T., Berends, C. J., Born, A., Box, J. E., Delhasse, A., Fujita, K., Gierz, P., Goelzer, H., Hanna, E., Hashimoto, A., Huybrechts, P., Kapsch, M.-L., King, M. D., Kittel, C., Lang, C., Langen, P. L., Lenaerts, J. T. M., Liston, G. E., Lohmann, G., Mernild, S. H., Mikolajewicz, U., Modali, K., Mottram, R. H., Niwano, M., Noël, B., Ryan, J. C., Smith, A., Streffing, J., Tedesco, M., van de Berg, W. J., van den Broeke, M., van de Wal, R. S. W., van Kampenhout, L., Wilton, D., Wouters, B., Ziemen, F., and Zolles, T.: GrSMBMIP: intercomparison of the modelled 1980–2012 surface mass balance over the Greenland Ice Sheet, The Cryosphere, 14, 3935–3958,, 2020. 

Forster, R. R., Box, J. E., van den Broeke, M. R., Miège, C., Burgess, E. W., van Angelen, J. H., Lenaerts, J. T. M., Koenig, L. S., Paden, J., Lewis, C., Gogineni, S. P., Leuschen, C., and McConnell, J. R.: Extensive liquid melt water storage in firn within the Greenland ice sheet, Nat. Geosci., 7, 95–98,, 2013. 

Giese, A. L. and Hawley, R. L.: Reconstructing thermal properties of firn at Summit, Greenland, from a temperature profile time series, J. Glaciol., 61, 503–510,, 2015. 

Griffiths, T. M.: Glaciological investigations in the TUTO area of Greenland., U. S. Army Snow Ice and Permafrost Research Establishment, Corps of Engineers, Report 47, 62 pp., 1960. 

Hanna, E., Cappelen, J., Fettweis, X., Mernild, S. H., Mote, T. L., Mottram, R., Steffen, K., Ballinger, T. J., and Hall, R. J.: Greenland surface air temperature changes from 1981 to 2019 and implications for ice-sheet melt and mass-balance change, Int. J. Climatol., 41, E1336–E1352,, 2021. 

Hansen, N., Langen, P. L., Boberg, F., Forsberg, R., Simonsen, S. B., Thejll, P., Vandecrux, B., and Mottram, R.: Downscaled surface mass balance in Antarctica: impacts of subsurface processes and large-scale atmospheric circulation, The Cryosphere, 15, 4315–4333,, 2021. 

Harrington, J. A., Humphrey, N. F., and Harper, J. T.: Temperature distribution and thermal anomalies along a flowline of the Greenland ice sheet, Ann. Glaciol., 56, 98–104,, 2015. 

Heilig, A., Eisen, O., MacFerrin, M., Tedesco, M., and Fettweis, X.: Seasonal monitoring of melt and accumulation within the deep percolation zone of the Greenland Ice Sheet and comparison with simulations of regional climate modeling, The Cryosphere, 12, 1851–1866,, 2018. 

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049,, 2020. 

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., and Thépaut, J.-N.: ERA5 monthly averaged data on single levels from 1940 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set],, 2023. 

Heuberger, J.-C.: Volume 1: Forages sur l'inlandsis, in: Expéditions Polaires Françaises: Missions Paul-Émile Victor, Glaciologie Groenland, Hermann and Cie, Paris, 1954. 

Hills, B. H., Harper, J. T., Meierbachtol, T. W., Johnson, J. V., Humphrey, N. F., and Wright, P. J.: Processes influencing heat transfer in the near-surface ice of Greenland's ablation zone, The Cryosphere, 12, 3215–3227,, 2018. 

Hirashima, H., Yamaguchi, S., Sato, A., and Lehning, M.: Numerical modeling of liquid water movement through layered snow based on new measurements of the water retention curve, Cold Reg. Sci. Technol., 64, 94–103,, 2010. 

Hofer, S., Tedstone, A. J., Fettweis, X., and Bamber, J. L.: Decreasing cloud cover drives the recent mass loss on the Greenland ice sheet, Sci. Adv., 3, e1700584,, 2017. 

Horlings, A. N., Christianson, K., and Miège, C.: Expansion of firn aquifers in southeast Greenland, J. Geophys. Res.-Earth Surf., 127, e2022JF006753,, 2022. 

How, P., Abermann, J., Ahlstrøm, A. P., Andersen, S. B., Box, J. E., Citterio, M., Colgan, W. T., Fausto, R. S., Karlsson, N. B., Jakobsen, J., Langley, K., Larsen, S. H., Mankoff, K. D., Pedersen, A. Ø., Rutishauser, A., Shield, C. L., Solgaard, A. M., van As, D., Vandecrux, B., and Wright, P. J.: PROMICE and GC-Net automated weather station data in Greenland, GEUS Dataverse [data set],, 2022. 

Howat, I. M., Negrete, A., and Smith, B. E.: The Greenland Ice Mapping Project (GIMP) land classification and surface elevation data sets, The Cryosphere, 8, 1509–1518,, 2014. 

Humphrey, N., Harper, J., and Meierbachtol, T.: Physical limits to meltwater penetration in firn, J. Glaciol., 67, 952–960,, 2021. 

Humphrey, N. F., Harper, J. T., and Pfeffer, W. T.: Thermal tracking of meltwater retention in Greenland's accumulation area, J. Geophys. Res.-Earth Surf., 117, F01010,, 2012. 

Kingma, D. P. and Ba, J.: Adam: A method for stochastic optimization, arXiv preprint arXiv: 1412.6980, 2014. 

Kjær, H. A., Zens, P., Edwards, R., Olesen, M., Mottram, R., Lewis, G., Terkelsen Holme, C., Black, S., Holst Lund, K., Schmidt, M., Dahl-Jensen, D., Vinther, B., Svensson, A., Karlsson, N., Box, J. E., Kipfstuhl, S., and Vallelonga, P.: Recent North Greenland temperature warming and accumulation, The Cryosphere Discuss. [preprint],, 2021. 

Koch, J. P. and Wegener, A.: Wissenschaftliche Ergebnisse Der Dänischen Expedition Nach Dronning Louises-Land Und Quer über Das Inlandeis Von Nordgrönland 1912–13 Unter Leitung Von Hauptmann J. P. Koch, Meddelelser om Grönland, Bd. 75, edited by: Reitzel, C. A., Copenhagen, 1930. 

Kuipers Munneke, P., van den Broeke, M. R., Reijmer, C. H., Helsen, M. M., Boot, W., Schneebeli, M., and Steffen, K.: The role of radiation penetration in the energy budget of the snowpack at Summit, Greenland, The Cryosphere, 3, 155–165,, 2009. 

Kuipers Munneke, P., Ligtenberg, S. R. M., van den Broeke, M. R., van Angelen, J. H., and Forster, R. R.: Explaining the presence of perennial liquid water bodies in the firn of the Greenland ice sheet, Geophys. Res. Lett., 41, 476–483,, 2014. 

Lance, G. N. and Williams, W. T.: Mixed-Data Classificatory Programs I – Agglomerative Systems, Australian Comput. J., 1, 15–20, 1967. 

Langen, P. L., Fausto, R. S., Vandecrux, B., Mottram, R. H., and Box, J. E.: Liquid Water Flow and Retention on the Greenland Ice Sheet in the Regional Climate Model HIRHAM5: Local and Large-Scale Impacts, Front. Earth Sci., 4, 110,, 2017. 

Langway, C. C.: Accumulation and Temperature on The Inland Ice of North Greenland, 1959, J. Glaciol., 3, 1017–1044,, 1961. 

Laternser, M.: Firn temperature and snow pit studies on the EGIG eastern traverse of central Greenland, 1992. Zürich, Eidgenössische Technische Hochschute, Versuchsanstalt Kur Wasserbau, Hydrologie und Glaziologie, Arbeitsheft 15, 1994. 

Law, R., Christoffersen, P., Hubbard, B., Doyle, S. H., Chudley, T. R., Schoonman, C. M., Bougamont, M., des Tombe, B., Schilperoort, B., Kechavarzi, C., Booth, A., and Young, T. J.: Thermodynamics of a fast-moving Greenlandic outlet glacier revealed by fiber-optic distributed temperature sensing, Sci. Adv., 7, 7136–7150,, 2021. 

Ligtenberg, S. R. M., Helsen, M. M., and van den Broeke, M. R.: An improved semi-empirical model for the densification of Antarctic firn, The Cryosphere, 5, 809–819,, 2011. 

Ligtenberg, S. R. M., Kuipers Munneke, P., Noël, B. P. Y., and van den Broeke, M. R.: Brief communication: Improved simulation of the present-day Greenland firn layer (1960–2016), The Cryosphere, 12, 1643–1649,, 2018. 

Liu, Q., Niu, J., Lu, P., Dong, F., Zhou, F., Meng, X., Xu, W., Li, S., and Hu, B. X.: Interannual and seasonal variations of permafrost thaw depth on the Qinghai-Tibetan plateau: A comparative study using long short-term memory, convolutional neural networks, and random forest, Sci. Total Environ., 838, 155886,, 2022. 

Loeb, N. A., Crawford, A., Stroeve, J. C., and Hanesiak, J.: Extreme precipitation in the Eastern Canadian Arctic and Greenland: An evaluation of atmospheric reanalyses, Front. Environ. Sci., 10, 866929,, 2022. 

Løkkegaard, A., Mankoff, K. D., Zdanowicz, C., Clow, G. D., Lüthi, M. P., Doyle, S. H., Thomsen, H. H., Fisher, D., Harper, J., Aschwanden, A., Vinther, B. M., Dahl-Jensen, D., Zekollari, H., Meierbachtol, T., McDowell, I., Humphrey, N., Solgaard, A., Karlsson, N. B., Khan, S. A., Hills, B., Law, R., Hubbard, B., Christoffersen, P., Jacquemart, M., Seguinot, J., Fausto, R. S., and Colgan, W. T.: Greenland and Canadian Arctic ice temperature profiles database, The Cryosphere, 17, 3829–3845,, 2023. 

Lorentzen, M., Bredesen, K., Mosegaard, K., and Nielsen, L.: Estimation of shear sonic logs in the heterogeneous and fractured Lower Cretaceous of the Danish North Sea using supervised learning, Geophys. Prospect., 70, 1410–1431,, 2022. 

Lundin, J. M. D., Stevens, C. M., Arthern, R., Buizert, C., Orsi, A., Ligtenberg, S. R. M., Simonsen, S. B., Cummings, E., Essery, R., Leahy, W., Harris, P., Helsen, M. M., and Waddington, E. D.: Firn Model Intercomparison Experiment (FirnMICE), J. Glaciol., 63, 401–422,, 2017. 

Lüthi, M. P., Ryser, C., Andrews, L. C., Catania, G. A., Funk, M., Hawley, R. L., Hoffman, M. J., and Neumann, T. A.: Heat sources within the Greenland Ice Sheet: dissipation, temperate paleo-firn and cryo-hydrologic warming, The Cryosphere, 9, 245–253,, 2015. 

MacFerrin, M., Machguth, H., As, D.V., Charalampidis, C., Stevens, C.M., Heilig, A., Vandecrux, B., Langen, P. L., Mottram, R., Fettweis, X., and van den Broeke, M. R.: Rapid expansion of Greenland's low-permeability ice slabs, Nature, 573, 403–407,, 2019. 

MacFerrin, M., Stevens, C. M., and Vandecrux, B.: The Greenland Firn Compaction Verification and Reconnaissance (FirnCover) Dataset, 2013–2019, Arctic Data Center [data set],, 2021. 

MacFerrin, M. J., Stevens, C. M., Vandecrux, B., Waddington, E. D., and Abdalati, W.: The Greenland Firn Compaction Verification and Reconnaissance (FirnCover) dataset, 2013–2019, Earth Syst. Sci. Data, 14, 955–971,, 2022. 

Machguth, H., MacFerrin, M., van As, D., Box, J. E., Charalampidis, C., Colgan, W., Fausto, R. S., Meijer, H. A., Mosley-Thompson, E., and van de Wal, R. S.: Greenland meltwater storage in firn limited by near-surface ice formation, Nat. Clim. Change, 6, 390–393,, 2016. 

Mankoff, K., Løkkegaard, A., Colgan, W., Thomsen, H., Clow, G., Fisher, D., Zdanowicz, C., Lüthi, M.P., Vinther, B., MacGregor, J.A., McDowell, I., Zekollari, H., Meierbachtol, T., Doyle, S., Law, R., Hills, B., Harper, J., Humphrey, N., Hubbard, B., Christoffersen, P., and Jacquemart, M.: Greenland deep ice temperature database, V1, GEUS Dataverse [data set],, 2022. 

Marshall, S. J.: Regime Shifts in Glacier and Ice Sheet Response to Climate Change: Examples From the Northern Hemisphere, Front. Clim., 3, 702585,, 2021. 

Matoba, S., Motoyama, H., Fujita, K., Yamasaki, T., Minowa, M., Onuma, Y., Komuro, Y., Aoki, T., Yamaguchi, S., Sugiyama, S., and Enomoto, H.: Glaciological and meteorological observations at the SIGMA-D site, northwestern Greenland ice sheet, Bull. Glaciol. Res., 33, 7–14,, 2015. 

McGrath, D., Colgan, W., Bayou, N., Muto, A., and Steffen, K.: Recent warming at Summit, Greenland: Global context and implications, Geophys. Res. Lett., 40, 2091–2096,, 2013. 

Meier, M. F., Conel, J. E., Hoerni, J. A., Melbourne, W. G., and Pings, C. J.: Preliminary Study of Crevasse Formation. Blue Ice Valley, Greenland, 1955, U.S. Army Snow Ice and Permafrost Research Establishment, Corps of Engineers, Report 38, 80 pp., 1957. 

Meyer, H. and Pebesma, E.: Predicting into unknown space? Estimating the area of applicability of spatial prediction models, Meth. Ecol. Evol., 12, 1620–1633,, 2021. 

Miller, N. B., Shupe, M. D., Cox, C. J., Noone, D., Persson, P. O. G., and Steffen, K.: Surface energy budget responses to radiative forcing at Summit, Greenland, The Cryosphere, 11, 497–516,, 2017. 

Miller, O., Solomon, D. K., Miège, C., Koenig, L., Forster, R., Schmerr, N., Ligtenberg, S. R. M., Legchenko, A., Voss, C. I., Montgomery, L., and McConnell, J. R.: Hydrology of a Perennial Firn Aquifer in Southeast Greenland: An Overview Driven by Field Data, Water Resour. Res., 56, e2019WR026348,, 2020. 

Mock, S. J.: Glaciological studies in the vicinity of Camp Century, Greenland, U.S. Army Cold Region Research and Engineering Report no. 157, (last access: 26 August 2022), 1965. 

Mock, S. J. and Ragle, R. H.: Elevations on the ice sheet of southern Greenland, U.S. Army Cold Region Research and Engineering Report no. 124, (last access: 26 August 2022), 1963. 

Mock, S. J. and Weeks, W. F.: The distribution of ten-meter snow temperatures on the Greenland ice sheet, U.S. Army Cold Region Research and Engineering Report no. 157, (last access: 26 August 2022), 1965. 

Muñoz Sabater, J.: ERA5-Land monthly averaged data from 1950 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS) [data set],, 2019. 

Nobles, L. H.: Glaciological investigations, Nunatarssuaq ice ramp, Northwestern Greenland, Tech. Rep. 66, U.S. Army Snow, Ice and Permafrost Research Establishment, Corps of Engineers, (last access: 26 August 2022), 1960. 

Noël, B., van de Berg, W. J., Lhermitte, S., and van den Broeke, M. R.: Rapid ablation zone expansion amplifies north Greenland mass loss, Sci. Adv., 5, 2–11,, 2019. 

Nolin, A. and Stroeve, J.: The changing albedo of the Greenland ice sheet: Implications for climate modeling, Ann. Glaciol., 25, 51–57,, 1997. 

Ohmura, A., Steffen, K., Blatter, H., Greuell, W., Rotach, M., Stober, M., Konzelmann, T., Forrer, J., Abe-Ouchi, A., Steiger, D., and Niederbaumer, G.: Energy and mass balance during the melt season at the equilibrium line altitude, Paakitsoq, Greenland ice sheet: Progress report, 2, 1992. 

Pfeffer, W. T., Meier, M. F., and Illangasekare, T. H.: Retention of Greenland runoff by refreezing: implications for projected future sea level change, J. Geophys. Res., 96, 22117,, 1991. 

Phillips, T., Rajaram, H., and Steffen, K.: Cryo – hydrologic warming: A potential mechanism for rapid thermal response of ice sheets, Geophys. Res. Lett., 37, 1–5,, 2010. 

Phillips, T., Rajaram, H., Colgan, W., Steffen, K., and Abdalati, W.: Evaluation of cryo-hydrologic warming as an explanation for increased ice velocities in the wet snow zone, Sermeq Avannarleq, West Greenland, J. Geophys. Res.-Earth Surf., 118, 1241–1256,, 2013. 

Poinar, K., Joughin, I., Lenaerts, J. T., and Van Den Broeke, M. R.: Englacial latent-heat transfer has limited influence on seaward ice flux in western Greenland, J. Glaciol., 63, 1–16,, 2017. 

Polashenski, C., Courville, Z., Benson, C., Wagner, A., Chen, J., Wong, G., Hawley, R., and Hall, D.: Observations of pronounced Greenland ice sheet firn warming and implications for runoff production, Geophys. Res. Lett., 41, 4238–4246,, 2014. 

Ragle, R. H. and Davis, T. C.: South Greenland traverses, J. Glaciol., 4, 129–131,, 1962. 

Rantanen, M., Karpechko, A. Y., Lipponen, A., Nordling, K., Hyvärinen, O., Ruosteenoja, K., Vihma, T., and Laaksonen, A.: The Arctic has warmed nearly four times faster than the globe since 1979, Nat. Commun. Earth Environ., 3, 168,, 2022. 

Rinker, J. and Mock, S.: Radar Ice Thickness Profiles: Northwest Greenland, CRREL Spec. Rep., 103, (last access: 26 August 2022), 1967. 

Rumelhart, D. E., Hinton, G. E., and Williams, R. J.: Learning representations by backpropagating errors, Nature, 323, 533–536,, 1986. 

Ryan, J. C., Smith, L. C., Van As, D., Cooley, S. W., Cooper, M. G., Pitcher, L. H., and Hubbard, A.: Greenland ice sheet surface melt amplified by snowline migration and bare ice exposure, Sci. Adv., 5, eaav3738,, 2019. 

Ryan, J. C., Smith, L. C., Wu, M., Cooley, S. W., Miège, C., Montgomery, L. N., Koenig, L. S., Fettweis, X., Noel, B. P. Y., and van den Broeke, M. R.: Evaluation of Cloudsat's cloud-profiling radar for mapping snowfall rates across the Greenland Ice Sheet, J. Geophys. Res.-Atmos., 125, e2019JD031411,, 2020. 

Ryan, J. C., Smith, L. C., Cooley, S. W., Pearson, B., Wever, N., Keenan, E., and Lenaerts, J. T. M.: Decreasing surface albedo signifies a growing importance of clouds for Greenland Ice Sheet meltwater production, Nat. Commun., 13, 4205,, 2022. 

Samimi, S., Marshall, S. J., Vandecrux, B., and MacFerrin, M.: Time-domain reflectometry measurements and modeling of firn meltwater infiltration at DYE-2, Greenland, J. Geophys. Res.-Earth Surf., 126, e2021JF006295,, 2021. 

Schwager, M.: Eisbohrkernuntersuchungen zur räumlichen und zeitlichen Variabilität von Temperatur und Niederschlagsrate im Spätholozän in Nordgrönland (Ice core analysis on the spatial and temporal variability of temperature and precipitation during the late Holocene in North Greenland), in: Berichte zur Polarforschung, Reports on Polar Research, Alfred-Wegener-Institut für Polar- und Meeresforschung, 362, 1–136,, 2000. 

Schytt, V.: Glaciological investigations in the Thule Ramp area, U. S. Army Snow Ice and Permafrost Research Establishment, Corps of Engineers, Report 28, 88 pp., (last access: 7 February 2024), 1955. 

Smeets, P. C., Kuipers Munneke, J. P., P., van As, D., van den Broeke, M. R., Boot, W., Oerlemans, H., Snellen, H., Reijmer, C. H., and van de Wal, R. S. W.: The K-transect in west Greenland: Automatic weather station data (1993–2016), Arct. Antarct. Alp. Res., 50, S100002,, 2018. 

Sorge, E.: Glaziologische Untersuchungen in Eismitte, in K. Wegener: Wissenschaftliche Ergebnisse der Deutschen Grönlandexpedition Alfred Wegener 1929 und 1930/1931, F. A. Brockhaus, 3, 1935. 

Stauffer, B. and Oeschger, H.: Temperaturprofile in Bohrlöchern am Rande des grönländischen Inlandeises, Mitteilungen der Versuchsanstalt für Wasserbau, Hydrologie und Glaziologie an der Eidgenössischen Technischen Hochschule (Zürich), 41, 301–313, 1979. 

Steffen, K. and Box, J.: Surface climatology of the Greenland ice sheet: Greenland Climate Network 1995–1999, J. Geophys. Res., 106, 33951–33964,, 2001. 

Steffen, K., Box, J., and Abdalati, W.: Greenland Climate Network: GC-Net, in: Glaciers, Ice Sheets and Volcanoes: A Tribute to Mark F. Meier, edited by: Colbeck, C., Special Report, Cold Regions Research and Engineering Laboratory, 1996. 

Steffen, K., Vandecrux, B., Houtz, D., Abdalati, W., Bayou, N., Box, J., Colgan, L., Espona Pernas, L., Griessinger, N., Haas-Artho, D., Heilig, A., Hubert, A., Iosifescu Enescu, I., Johnson-Amin, N., Karlsson, N. B., Kurup Buchholz, R., McGrath, D., Cullen, N. J., Naderpour, R., Molotch, N. P., Pederson, A. Ø., Perren, B., Philipps, T., Plattner, G. K., Proksch, M., Revheim, M. K., Særrelse, M., Schneebli, M., Sampson, K., Starkweather, S., Steffen, S., Stroeve, J., Watler, B., Winton, Ø. A., Zwally, J., and Ahlstrøm, A.: GC-Net Level 1 automated weather station data, V3, GEUS Dataverse [data set],, 2023. 

Steger, C. R., Reijmer, C. H., and van den Broeke, M. R.: The modelled liquid water balance of the Greenland Ice Sheet, The Cryosphere, 11, 2507–2526,, 2017. 

Steiner, D., Walter, A., and Zumbühl, H.: The application of a non-linear back-propagation neural network to study the mass balance of Grosse Aletschgletscher, Switzerland, J. Glaciol., 51, 313–323,, 2005. 

Stibal, M., Box, J. E., Cameron, K. A., Langen, P. L., Yallop, M. L., Mottram, R. H., Khan, A. L., Molotch, N. P., Chrismas, N. A. M., Calì Quaglia, F., Remias, D., Smeets, C. J. P. P., van den Broeke, M. R., Ryan, J. C., Hubbard, A., Tranter, M., van As, D., and Ahlstrøm, A. P.: Algae Drive Enhanced Darkening of Bare Ice on the Greenland ice sheet, Geophys. Res. Lett., 44, 11463–11471,, 2017. 

Tedstone, A. J. and Machguth, H.: Increasing surface runoff from Greenland's firn areas, Nat. Clim. Chang., 12, 672–676,, 2022. 

Thomsen, H., Olesen, O., Braithwaite, R., and Bøggild, C.: Ice drilling and mass balance at Pâkitsoq, Jakobshavn, central West Greenland, Rapp. Grønlands Geol. Undersøgelse, 152, 80–84,, 1991. 

Trusel, L. D., Das, S. B., Osman, M. B., Evans, M. J., Smith, B. E., Fettweis, X., McConnell, J. R., Noël, B. P., and van den Broeke, M. R.: Nonlinear rise in Greenland runoff in response to post-industrial Arctic warming, Nature, 564, 104–108,, 2018. 

U.S. Army Transportation Board: Report of environmental operation: Lead Dog 1960, Final Report, Project Lead Dog, TCB- 60 - 023 - E0, 188 pp., (last access: 7 February 2024), 1960. 

van Dalum, C. T., van de Berg, W. J., and van den Broeke, M. R.: Impact of updated radiative transfer scheme in snow and ice in RACMO2.3p3 on the surface mass and energy budget of the Greenland ice sheet, The Cryosphere, 15, 1823–1844,, 2021. 

Vandecrux, B.: Greenland ice sheet 10 m subsurface temperature compilation 1912–2022, V1, GEUS Dataverse [data set],, 2023a. 

Vandecrux, B.: Greenland ice sheet 10 m ice and firn temperature reconstruction, 1950–2022, V2, GEUS Dataverse [data set],, 2023b. 

Vandecrux, B.: Greenland ice sheet subsurface temperature compilation scripts (v1.0), Zenodo [code],, 2023c. 

Vandecrux, B., MacFerrin, M., Machguth, H., Colgan, W. T., van As, D., Heilig, A., Stevens, C. M., Charalampidis, C., Fausto, R. S., Morris, E. M., Mosley-Thompson, E., Koenig, L., Montgomery, L. N., Miège, C., Simonsen, S. B., Ingeman-Nielsen, T., and Box, J. E.: Firn data compilation reveals widespread decrease of firn air content in western Greenland, The Cryosphere, 13, 845–859,, 2019. 

Vandecrux, B., Fausto, R. S., Van As, D., Colgan, W., Langen, P. L., Haubner, K., Ingeman-Nielsen, T., Heilig, A., Stevens, C. M., MacFerrin, M., Niwano, M., Steffen, K., and Box, J. E.: Firn cold content evolution at nine sites on the Greenland ice sheet between 1998 and 2017, J. Glaciol., 66, 591–602,, 2020a. 

Vandecrux, B., Mottram, R., Langen, P. L., Fausto, R. S., Olesen, M., Stevens, C. M., Verjans, V., Leeson, A., Ligtenberg, S., Kuipers Munneke, P., Marchenko, S., van Pelt, W., Meyer, C. R., Simonsen, S. B., Heilig, A., Samimi, S., Marshall, S., Machguth, H., MacFerrin, M., Niwano, M., Miller, O., Voss, C. I., and Box, J. E.: The firn meltwater Retention Model Intercomparison Project (RetMIP): evaluation of nine firn models at four weather station sites on the Greenland ice sheet, The Cryosphere, 14, 3785–3810,, 2020b. 

Vandecrux, B., Colgan, W., Solgaard, A. M., Steffensen, J. P., and Karlsson, N. B.: Firn evolution at Camp Century, Greenland: 1966–2100, Front. Earth Sci., 9,, 2021. 

Vandecrux, B., Amory, C., Ahlstrøm, A. P., Akers, P. D., Albert, M., Alley, R. B., Arnaud, L., Bales, R., Benson, C., Box, J. E., Buizert, C., Charalampidis, C., Clerx, N., Covi, F., Denis, G., Dibb, J. E., Ding, M., Eisen, O., Fausto, R., Fernandoy, F., Freitag, J., Gerland, S., Harper, J., Hawley, R. L., Hock, R., How, P., Hubbard, B., Humphrey, N., Iizuka, Y., Isaksson, E., Kameda, T., Karlsson, N. B., Kawakami, K., Kjær, H. A., Kuipers Munneke, P., Lewis, G., MacFerrin, M., Machguth, H., Mankoff, K. D., McConnell, J. R., Medley, B., Morris, E., Mosley-Thompson, E., Mulvaney, R., Niwano, M., Osterberg, E., Otosaka, I., Picard, G., Polashenski, C., Rennermalm, A., Rutishauser, A., Simonsen, S. B., Smith, A., Solgaard, A., Spencer, M., Steen-Larsen, H. C., Stevens, C. M., Sugiyama, S., Tedesco, M., Thompson-Munson, M., Tsutaki, S., van As, D., Van den Broeke, M. R., Wilhelms, F., Xiao, J., and Xiao, C.: The SUMup collaborative database: Surface mass balance, subsurface temperature and density measurements from the Greenland and Antarctic ice sheets (1912–2023), Arctic Data Center [data set],, 2023a. 

Vandecrux, B., Box, J. E., Ahlstrøm, A. P., Andersen, S. B., Bayou, N., Colgan, W. T., Cullen, N. J., Fausto, R. S., Haas-Artho, D., Heilig, A., Houtz, D. A., How, P., Iosifescu Enescu, I., Karlsson, N. B., Kurup Buchholz, R., Mankoff, K. D., McGrath, D., Molotch, N. P., Perren, B., Revheim, M. K., Rutishauser, A., Sampson, K., Schneebeli, M., Starkweather, S., Steffen, S., Weber, J., Wright, P. J., Zwally, H. J., and Steffen, K.: The historical Greenland Climate Network (GC-Net) curated and augmented level-1 dataset, Earth Syst. Sci. Data, 15, 5467–5489,, 2023b. 

van den Broeke, M., Smeets, P., Ettema, J., van der Veen, C., van de Wal, R., and Oerlemans, J.: Partitioning of melt energy and meltwater fluxes in the ablation zone of the west Greenland ice sheet, The Cryosphere, 2, 179–189,, 2008. 

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791,, 2012. 

Wang, W., Zender, C. S., van As, D., Fausto, R. S., and Laffin, M. K.: Greenland surface melt dominated by solar and sensible heating, Geophys. Res. Lett., 48, e2020GL090653,, 2021. 

Weertman, J.: Comparison between measured and theoretical temperature profiles of the Camp Century, Greenland, Borehole, J. Geophys. Res., 73, 2691–2700,, 1968. 

Wegener, A.: Deutsche Grönland- Expedition Alfred Wegener: ausgeführt unter der Leitung von Alfred Wegener und Kurt Wegener 1929 und 1930/1931: wissenschaftliche Ergebnisse, 1940. 

Xiong, Z., Cui, Y., Liu, Z., Zhao, Y., Hu, M., and Hu, J.: Evaluating explorative prediction power of machine learning algorithms for materials discovery using k-fold forward cross-validation, Comput. Mater. Sci., 171, 109203,, 2020. 

Xu, W., Balaguru, K., August, A., Lalo, N., Hodas, N., DeMaria, M., and Judi, D.: Deep Learning Experiments for Tropical Cyclone Intensity Forecasts, Weather Forecast., 36, 1453–1470,, 2021. 

Yamaguchi, S., Matoba, S., Yamazaki, T., Tsushima, A., Niwano, M., Tanikawa, T., and Aoki, T.: Glaciological observations in 2012 and 2013 at SIGMA-A site, Northwest Greenland, Bull. Glaciol. Res., 32, 95–105,, 2014. 

Yen, Y.-C.: Review of thermal properties of snow, ice and sea ice, CRREL Report 81-10, 1–27, (last access: 2 November 2020), 1981. 

Zhang, Q., Huai, B., van den Broeke, M. R., Cappelen, J., Ding, M., Wang, Y., and Sun, W.: Temporal and spatial variability in contemporary Greenland warming (1958–2020), J. Climate, 35, 2755–2767,, 2022. 

Zwally, H. J., Giovinetto, M. B., Beckley, M. A., and Saba, J. L.: Antarctic and Greenland Drainage Systems, GSFC Cryospheric Sciences Laboratory, (last access: 7 February 2024), 2012. 

Short summary
How fast is the Greenland ice sheet warming? In this study, we compiled 4500+ temperature measurements at 10 m below the ice sheet surface (T10m) from 1912 to 2022. We trained a machine learning model on these data and reconstructed T10m for the ice sheet during 1950–2022. After a slight cooling during 1950–1985, the ice sheet warmed at a rate of 0.7 °C per decade until 2022. Climate models showed mixed results compared to our observations and underestimated the warming in key regions.