# Estimating the upper depth of subsurface water on the Greenland Ice Sheet using multi-frequency passive microwave remote sensing, radiative transfer modeling, and machine learning

Baptiste Vandcrux, Ghislain Picard, Pierre Zeiger, Marion Leduc-Leballeur, Andreas Colliander, Alamgir Hossan, Andreas Ahlstrøm

<figure><img src="image_1.png" /></figure>

## To cite this version:

Baptiste Vandcrux, Ghislain Picard, Pierre Zeiger, Marion Leduc-Leballeur, Andreas Colliander, et al.. Estimating the upper depth of subsurface water on the Greenland Ice Sheet using multi-frequency passive microwave remote sensing, radiative transfer modeling, and machine learning. Remote Sensing of Environment, 2026, 334, pp.115197. 10.1016/j.rse.2025.115197. hal-05464656

HAL Id: hal-05464656

https://hal.science/hal-05464656v1

Submitted on 19 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# Estimating the upper depth of subsurface water on the Greenland Ice Sheet using multi-frequency passive microwave remote sensing, radiative transfer modeling, and machine learning

Baptiste Vandecrux<sup>a,*</sup>, Ghislain Picard<sup>b</sup>, Pierre Zeiger<sup>c</sup>, Marion Leduc-Leballeur<sup>d</sup>, Andreas Colliander<sup>e,f</sup>, Alangir Hossan<sup>e</sup>, Andreas Ahlstrøm<sup>a</sup>

## ARTICLE INFO

Editor: Dr. Menghua Wang

**Keywords:**

Meltwater

Greenland

Ice sheet

Snow

Firn

Refreezing

## ABSTRACT

As the Arctic warms, surface melt extends into the Greenland Ice Sheet's accumulation zone, where much of the water infiltrates into the snowpack. This makes monitoring the subsurface water depth and spatial extent important for accurate ice sheet runoff estimations. Subsurface water can be detected using remotely sensed microwave brightness temperatures ($T_B$). We use vertically polarized $T_B$ at 1.4 GHz from Soil Moisture and Ocean Salinity satellite (SMOS) and at 6.9, 10.7, and 18.7 GHz from the Advanced Microwave Scanning Radiometers (AMSR-E/2) to estimate the upper depth of liquid water (UDLW) on the ice sheet accumulation area. We build a catalogue of simulated UDLW and $T_B$: realistic UDLW are modeled by the Geological Survey of Denmark and Greenland (GEUS) snow model, forced by the Copernicus Arctic Regional Reanalysis (CARRA), and the corresponding $T_B$ are calculated by the Snow Microwave Radiative Transfer (SMRT) model at 19 sites. We train on this catalogue an ensemble of cross-validated Random Forest (RF) models to predict UDLW and its uncertainty from $T_B$ at four frequencies. On hold-out modeled data and for water within 5 m of the surface, the RF ensemble achieves a median RMSE of 0.68 m and mean error of -0.09 m. Our retrieval, when applied to observed $T_B$, matches within 2 m UDLW inferred from subsurface temperature profiles down to 4-6 m depth. Performances decrease beyond 5 m depth and for low liquid water amounts. Our retrieval produces daily UDLW maps over the ice sheet's accumulation area during 2010-2023 which reveal the seasonal evolution of UDLW, deliver the first quantitative estimates of subsurface liquid water depth on the ice sheet and offer new insights into meltwater infiltration and storage processes.

## 1. Introduction

Surface melt and runoff are responsible for about half of the Greenland Ice Sheet mass loss and contribution to sea-level rise (Otosaka et al., 2023). About 78 % of the Greenland Ice Sheet's runoff currently originates from the bare ice area (Steger et al., 2017), where there is enough energy input to the surface to melt the snow that accumulated over the previous winter and some of the underlying glacial ice. Above that region lies the accumulation area, covering 88 % of the ice sheet (Box et al., 2022), where melt does not exceed the amount of snow

accumulated during winter. The accumulation area has been subject to an increase in annual surface melt over recent decades (e.g., van den Broeke et al., 2016), which, when infiltrating into the snow and underlying firn, has either refrozen into multi-meter thick ice layers called ice slabs (Machguth et al., 2016; MacFerrin et al., 2019), superimposed ice (Tedstone et al., 2025) or contributed to the expansion of firn aquifers (Horlings et al., 2022). These changes in the infiltration of meltwater into the snow have also shifted the runoff zone – the area where surface meltwater initiates a downslope movement and eventually reaches the ocean – to higher elevations on the ice sheet (Tedstone

<footnote><sup>a</sup> Geological Survey of Denmark and Greenland, Copenhagen, Denmark</footnote>

<footnote><sup>b</sup> Univ. Grenoble Alpes, CNRS, Institut des Géosciences de l'Environnement (IGE), UMR 5001, Grenoble, France</footnote>

<footnote><sup>c</sup> CNES, Institut des Géosciences de l'Environnement, UMR 5001, Grenoble, France</footnote>

<footnote><sup>d</sup> Institute of Applied Physics "Nello Carrara", National Research Council, 50019 Sesto Fiorentino, Italy</footnote>

<footnote><sup>e</sup> Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, USA</footnote>

<footnote><sup>f</sup> Finnish Meteorological Institute, Helsinki, Finland</footnote>

<footnote>* Corresponding author at: Geological Survey of Denmark and Greenland, Øster Voldgade 10, 1350 Copenhagen, Denmark. E-mail address: bav@geus.dk (B. Vandecrux).</footnote>
and Machguth, 2022). While Tedstone and Machguth (2022) focused on visible surface features of runoff, a later study expended these results uphill of the visible runoff line and showed that meltwater infiltrates below the snow surface and refreezes either in the snow or as superimposed ice (Tedstone et al., 2025). These efforts were based on field observations from the so-called ice slab area in southwest Greenland and little is known about meltwater infiltrations in other regions of the ice sheet except for punctual measurements (e.g. Humphrey et al., 2012; Miller et al., 2020). A better understanding of the fate of the meltwater in the ice sheet accumulation area is needed for accurately predicting the ice sheet response to a changing climate (van den Broeke et al., 2023). Yet, to this day, there is no large-scale, quantitative observation of meltwater below the snow surface in the accumulation area of the ice sheet.

Passive microwave remote sensing, through the measurement of the brightness temperature ($T_B$) has been used for decades to evaluate the spatial and temporal extent of surface melt on ice sheets (e.g., Mote et al., 1993; Leduc-Leballeur et al., 2020; Colosio et al., 2021). The microwave thermal emission of the snow depends on its physical temperature, density, grain size, and ice content profiles (e.g., Zwally, 1977), as well as on the amount and depth of liquid water (e.g., Picard et al., 2022a). This sensitivity of $T_B$ to liquid water has been used to map the presence or absence of liquid water at the surface (e.g., Mote et al., 1993; Colosio et al., 2021) and, at lower frequencies, to quantify the total liquid water content of the snow (Houtz et al., 2021; Mousavi et al., 2021; Hossan et al., 2025). $T_B$ variations are representative of subsurface changes occurring down to a certain depth, increasing with decreasing frequency (Picard et al., 2022a; Colliander et al., 2022, 2023). Consequently, multi-frequency $T_B$ observations have been used to study subsurface liquid water at depths ranging from tens of centimeters for higher frequencies such as 18.7 or 36.5 GHz, to several meters for low frequencies such as 1.4 GHz, either through field observation at a single site in Greenland (Colliander et al., 2022) or through the sensitivity analysis of a microwave emission model in Antarctica (Picard et al., 2022a). Based on this qualitative understanding of the depth range to which each microwave frequency is sensitive, multi-frequency melt classifications have been developed for Greenland and Antarctica (Colliander et al., 2023; Leduc-Leballeur et al., 2025). These classifications indicate whether water is present at the surface based on higher frequencies such as 18.7 and 36.5 GHz or buried under a dry layer based on low frequencies such as 1.4 and 6.9 GHz, but do not attempt to quantify the depth at which this liquid water is buried. Due to the high absorptivity of wet snow in the microwave domain, the emission of a wet layer in the snowpack increasingly blocks most of the emission originating from deeper layers as the wetness increases. This means that multi-frequency microwave observations can detect the upper depth of liquid water (UDLW), but are unlikely to resolve its bottom depth or total thickness. Recently, this characteristic was used by Shang et al. (2025) to extrapolate the UDLW beyond airborne radar data acquired between 2011 and 2017 using spaceborne multi-frequency microwave observations. Developing this approach into an operational mapping of UDLW would be highly valuable for 1) providing information on the meltwater infiltration depth; 2) studying the location and temporal evolution of areas where subsurface meltwater is retained before it either refreezes or runs off; 3) preventing the misclassification of a deep wet layer as active surface melt. The relationship between $T_B$ observations at multiple frequencies and UDLW is highly non-linear and high-dimensional. While traditional retrievals of liquid water on the ice sheet have relied on model inversion (Houtz et al., 2021) or lookup tables (Hossan et al., 2025), machine learning provides a flexible alternative, as demonstrated by studies retrieving, e.g., melt intensity (Zheng et al., 2022) and snow characteristics (Chen et al., 2025) from satellite measurements.

The aim of this study is to provide a quantitative estimate of UDLW over the Greenland Ice Sheet using multi-frequency microwave observations and machine learning, thereby improving the description of

subsurface water, which until now has only been available in categorical form (Colliander et al., 2023; Leduc-Leballeur et al., 2025). Understanding the link between multi-frequency microwave observations and the profile of liquid water content in the snowpack is challenging when no large-scale observation of UDLW is available. To circumvent this lack of observations, we build a catalogue of simulated snow profiles containing liquid water from which both UDLW and multi-frequency $T_B$ are calculated. This catalogue is built as follows. At 19 well-documented sites, the multilayer snow and firn model developed at the Geological Survey of Denmark and Greenland (GEUS), henceforth called the GEUS snow model (Vandecrux et al., 2020a), driven by the Copernicus Arctic Regional Reanalysis (CARRA) dataset (Schyberg et al., 2020), is used to produce typical profiles of snow density, temperature and liquid water content, from which the UDLW can be calculated. These snow profiles are then used as inputs for the Snow Microwave Radiative Transfer (SMRT) model (Picard et al., 2018) to calculate the $T_B$ at 1.4, 6.9, 10.7, and 18.7 GHz. These simulations give a catalogue of UDLW and their corresponding $T_B$ on which we train an ensemble of cross-validated random forest (RF) models, which relate the catalogued multi-frequency $T_B$ to their corresponding UDLW. Last, we use this RF model ensemble on spaceborne $T_B$ observations at the same frequencies to provide an estimation of UDLW in the ice sheet's accumulation area for 2010–2023.

## 2. Data

### 2.1. Multi-frequency passive microwave observations

We use observations from the radiometers onboard three satellites. The Soil Moisture and Ocean Salinity (SMOS) satellite was launched by the European Space Agency (ESA) in 2009 and measures $T_B$ at 1.4 GHz with an average native spatial resolution of 43 km (Kerr et al., 2001). We here use the product from Zeiger et al. (2024), which harmonized the $T_B$ to a 40° incidence angle and enhanced the effective resolution to ~32 km. For higher frequencies, we use the two Advanced Microwave Scanning Radiometers: AMSR-E, onboard the NASA satellite Aqua, active during 2002–2011, and AMSR2, onboard the JAXA Global Change Observation Mission – Water (GCOM-W1) satellite, active since 2012. Both sensors, henceforth referred to as AMSR, measure $T_B$ at 6.9, 10.7, 18.7, 23.8, 36.5, and 89.0 GHz and at a 55° incidence angle. Here, we use the $T_B$ at 6.9, 10.7, and 18.7 GHz which have a spatial resolution of approximately 46, 32 and 18 km respectively. These frequencies are known to be sensitive to the water presence at various depths (Picard et al., 2022a; Colliander et al., 2022, 2023; Leduc-Leballeur et al., 2025) without being overly sensitive to surface properties as even higher frequencies can be.

The SMOS and AMSR products both provide horizontally and vertically polarized $T_B$, henceforth H-pol and V-pol $T_B$. The H-pol $T_B$ shows the largest increase between dry conditions and in presence of liquid water, but it is known to be sensitive to steep density gradients within the snow, called layering, which is variable in time and space (e.g., Picard et al., 2022a). Contrastingly, the V-pol $T_B$ shows slightly lower increase in presence of liquid water, but has the advantage of being less sensitive to the layering at the satellites' observation angles and more closely related to the snow grain size, temperature, and liquid water profiles (e.g., Picard et al., 2022a; Leduc-Leballeur et al., 2017, 2025). Consequently, we here only use V-pol $T_B$ to make our retrieval more robust when applied in regions with unknown layering. All V-pol $T_B$ grids are resampled on a common 12.5 km polar stereographic grid using nearest neighbor interpolation, except for the 18.7 GHz $T_B$ which is already distributed at this resolution.

Our multi-frequency approach leverages subtle changes in $T_B$ at different frequencies to extract information on the depth of the liquid water. SMOS and AMSR provide measurements over Greenland twice daily, but we find that ascending and descending measurements could in many places be more than 1 K different, even during polar night when
there should be no diurnal cycle. To minimize the $T_B$ variations that are not related to liquid water dynamics and maximize the synchronicity of measurements from different satellites, we only use the descending passes for AMSR collected between 5 and 9 UTC and ascending passes for SMOS collected between 5 and 8 UTC. Additionally, we find that the $T_B$ at 1.4 GHz shows a random noise of ~3 K in winter in the dry snow area, even though they should be representative of the very stable deep ice sheet structure (e.g., Leduc-Leballeur et al., 2017). This noise, likely coming from anthropic radiofrequency interferences and accentuated by the instrumental design and resolution enhancement algorithm, caused large variability in our retrieval output. To make our retrieval more stable over time, we therefore smooth all $T_B$ using a Gaussian mean on a seven-day running window. Additionally, since the microwave emission from the bare ice area is markedly different that of the accumulation area (Miller et al., 2022) and is challenging for retrievals of liquid water (Hossan et al., 2025), we here only use observed $T_B$ from the accumulation area as delineated by Vandecrux et al. (2019).

## 2.2. Single-frequency passive microwave surface wetness indicators

While our retrieval is designed to be most sensitive to the subsurface water, other algorithms have been developed to provide a binary indicator of surface wetness (e.g., Colosio et al., 2021; Zeiger et al., 2024). To ensure that our retrieval's output is restricted to locations and periods where liquid water is present at or below the surface, we use the snow wetness indicator derived from H-pol $T_B$ at 1.4 GHz from Zeiger et al. (2024): if no liquid water is detected at 1.4 GHz H-pol $T_B$, then we assume that water is either absent or beyond the detection capacity of our retrieval. To further evaluate our retrieval, we use the single-frequency snow wetness indicators derived through adaptive thresholding of V-pol $T_B$ at 6.9 and 18.7 GHz (Torinesi et al., 2003; Leduc-Leballeur et al., 2025).

## 2.3. The CARRA reanalysis

CARRA is the latest climate reanalysis dataset from the European Centre for Medium-Range Weather Forecasts (Schyberg et al., 2020). Focusing on a domain including part of the Canadian Arctic, Greenland, Iceland and Svalbard, it assimilates numerous observations including, in Greenland, remotely sensed surface albedo, and pressure, temperature, humidity, and wind speed measured by on-ice weather stations (Fausto et al., 2021; Vandecrux et al., 2023). With a spatial resolution of 2.5 km and a better description of cold surfaces, CARRA is significantly better compared to in situ observations than previous reanalyses (Køltzow et al., 2022).

## 2.4. Subsurface temperature measurements and site selection

Apart from punctual observation campaigns (e.g., Heilig et al., 2018; Samimi et al., 2020), there are virtually no direct and continuous observations of the liquid water content and its vertical distribution on the Greenland Ice Sheet. Nevertheless, subsurface temperature measured by thermistor strings can indicate when the snow is at the melting point and most likely contains water. Such instrumentation has been deployed broadly, because of its simplicity and cost-effectiveness, on the Greenland Ice Sheet (e.g., Humphrey et al., 2012; Charalampidis et al., 2016; Miège et al., 2016; Miller et al., 2020; Fausto et al., 2021; Steffen et al., 1996; Vandecrux et al., 2023; Saito et al., 2024) and has been recently compiled and harmonized in the SUMup collaborative database along with their measurement depth (Vandecrux et al., 2024a, 2024b). For the construction of our retrieval, we selected 19 sites, covering most of the Greenland Ice Sheet accumulation area, and where subsurface temperatures are available to assess the GEUS snow model runs and UDLW retrieval (Fig. 1). Among these sites, six are located close to the equilibrium line and receive abundant melt, between 0.3 and 0.6 m water equivalent (w.e.) every year according to the GEUS snow model

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Sites where the GEUS and SMRT models are run to build a training catalogue of $T_B$ and corresponding UDLW.</figcaption></figure>

(FA-13, FA-15-1, FA-15-2, H2, H3, KAN_U), ten sites are located in the percolation area where between 0.03 and 0.27 m w.e. of melt occurs every year (CP1, DYE-2, Saddle, South Dome, T1, T2, T3, T4, NASA-SE, SIGMA-A) and three sites are from the dry snow area where the annual surface melt is lower than 0.01 m w.e. (NASA-U, Humboldt, CEN).

Although thermistor strings are available at our 19 sites over various periods, their suitability for evaluating our retrieval varies. Shortly after installation, the freshly drilled borehole can allow meltwater to infiltrate along the instrument cable (Humphrey et al., 2012), which can distort the recorded temperature profiles. After a few years, the borehole can be considered sealed, but the instruments may then lie too deep to capture meltwater infiltration. Based on a careful screening of our subsurface temperature dataset, we identified three site-year combinations that provide the most reliable conditions for our evaluation: FA-13 in 2013, FA-15-1 in 2015, and CP1 in 2021.

Since our retrieval is trained on data from 19 sites, we also assess how it performs in regions without training data. For this purpose, we examine a location in the Qassimiut aquifer region (61.63120°N, -46.53871°E; Fig. 1), an area identified as a perennial firn aquifer by Miège et al. (2016) but without any in situ temperature or water-content observations. At this site, we compare our retrieved UDLW to liquid-water profiles simulated by the GEUS snow model forced with CARRA, allowing us to evaluate both the realism and the usefulness of the retrieval away from our training sites.

# 3. Methods

There are no large-scale in situ observations of UDLW available for comparison with satellite data. So, to investigate how $T_B$ varies across frequencies in response to changes in UDLW, we use a model suite that generates a catalogue of plausible snow profiles, each with associated UDLW values and corresponding multi-frequency $T_B$. This model chain is run at our 19 selected sites for 2010–2023 and provides the training dataset to our machine learning model which eventually estimates UDLW from a set of observed multi-frequency $T_B$.

## 3.1. The GEUS snow model

The GEUS snow model was originally developed as a combination of
the Surface Energy Budget (SEB) model from Van As et al. (2005), with the subsurface scheme of HIRHAM (Langen et al., 2017). It then received numerous updates and was tested extensively at various sites on the Greenland Ice Sheet (Vandecrux et al., 2018, 2020a, 2020b). Here, the GEUS model takes as inputs snowfall, rainfall, temperature, relative humidity, wind speed, net shortwave radiative flux, and downward longwave radiative flux from CARRA at three-hourly temporal resolution. The surface module then calculates the energy balance at the surface, and the surface melt, if any. Then the subsurface part of the model handles the diffusion of temperature into the snow, its densification, the growth of snow grain and the infiltration of meltwater through the 100 layers of the model. For each site of interest (Fig. 1), the GEUS model is spun up five times over the 1990–2001 climate, then run for 1990–2023. The UDLW for the period 2010–2023 is extracted from the model output for each time stamp as the depth of the uppermost wetted layer.

## 3.2. The SMRT model

To link the changes in $T_B$ to changes in the snow characteristics, including UDLW, we apply the SMRT model (Picard et al., 2018) on the profiles generated by the GEUS snow model. SMRT takes as input the physical and microstructural properties of a multilayer snowpack: each layer's type – snow or ice – grain size, density and temperature. SMRT then solves the radiative transfer equation for each layer: the amount of radiation scattered or absorbed by the media and reflected or transmitted by the interfaces. SMRT includes several formulations of radiative transfer in dense media, including the Improved Born Approximation (IBA, Mätzler, 1998), which we use here. To compute the scattering coefficient, IBA relies on an analytical form of the media autocorrelation function, which is linked to a user-defined microstructural snow model. Here, we use the sticky hard sphere model with a stickiness of 0.17 (Löwe and Picard, 2015 and references therein) with the effective permittivity mixing formula of Polder and Van Santeen (1946). Contributions from the atmosphere were not considered in the radiative transfer modeling as they were shown to be minor at our frequencies of interest (e.g., Kontu and Pulliainen, 2008).

## 3.3. Optimization of the coupling between the SMRT and GEUS models

When building machine learning models, the similarity of the model inputs in the training data ($T_B$ modeled by SMRT) and the target data (here $T_B$ observed by SMOS and AMSR) is crucial (e.g., Meyer and Pebesma, 2021). To ensure this, we implement the following coupling of the GEUS model and SMRT by optimizing two parameters each year to minimize the difference between modeled and observed $T_B$ during winter, thereby achieving more realistic simulated $T_B$ over the following melting season (Fig. 2a). Before the optimization, we resample the pure snow density, the temperature and the grain size from the 100 dynamical layers of the GEUS snow model to 50 fixed depth intervals using weighted means. The ice content from the GEUS snow model is resampled to the 50 fixed depth intervals using weighted sums.

Yet the ice content resampled in the 50 fixed layers does not specify the number or thickness of ice layers in each depth interval. Preliminary tests showed that for a given vertical profile of ice content, concentrating the ice into many thin pure ice layers, or into few thick ice layers gives different $T_B$ at 1.4 GHz due to reflections at snow-ice interfaces. The first optimization parameter is therefore the snowpack's total number of pure ice layers $N_{ice}$ that will be defined in the SMRT input. We develop a method that distributes the ice content of the 50 fixed depth intervals between a user-defined number $N_{ice}$ of pure ice layers, while preserving the total amount and the vertical distribution of ice in the profile (Fig. 2b). The method takes as a starting point the number of depth intervals containing ice. If $N_{ice}$ equals this number, then one pure ice layer is inserted in each ice-containing depth interval. Ice layers are centered in their respective depth interval and their thickness is defined by the ice content in the depth interval. If $N_{ice}$ is higher than the initial

number of ice-containing depth intervals, the method iteratively splits the ice content within each of the depth intervals into thinner ice layers until there are $N_{ice}$ ice layers in the 50 depth intervals (Fig. 2b). Ice layers are assumed equally spaced and of identical thickness within each of the depth interval. If $N_{ice}$ is lower than the initial number of depth intervals containing ice, we start with one ice layer for each ice-containing depth interval, and then iteratively aggregate the smallest ice layer to the second smallest ice layer until there are $N_{ice}$ ice layers left. The modification of $N_{ice}$ enables us to find the snowpack description that leads to the best match between observed winter $T_B$, defined as the median value of V-pol $T_B$ at 1.4 GHz for the first 60 days of the year and the modeled winter $T_B$. For each site and for each year, we take the median winter profiles of snow density, temperature, and ice content during the first 60 days of the year and iteratively find the optimal $N_{ice}$ value, between 1 and 200 (secant method, 10 iterations), that produces the smallest difference between SMRT's simulated $T_B$ and the observed winter $T_B$ at 1.4 GHz (Fig. 2a). The Supplementary Fig. 1a shows this adjustment process at the site DYE-2 in 2021.

The second adjustment adapts the snow grain size calculated by the GEUS model before it is used by SMRT (Fig. 2a). The GEUS snow model assumes that snow grains have a diameter of 0.1 mm when deposited at the surface and subsequently grow as a function of the snow temperature and wetness after Brun (1989). This formulation is known to be not directly comparable to the snow grain size needed by electromagnetic models and therefore needs adaptation (e.g., Brucker et al., 2011; Langlois et al., 2012; Picard et al., 2022b). We define a multiplicative correction factor called $C_{GS}$, which is to be applied to the grain size profile given by the GEUS model and resampled in the 50 fixed depth intervals. For each year and each site, the observed winter V-pol $T_B$ at 6.9, 10.7, and 18.7 GHz are calculated as the median value for the first 60 days of the year. Eventually, we adjust $C_{GS}$ iteratively (secant method, 10 iterations) to minimize the mean difference between the observed and modeled winter V-pol $T_B$ at 6.9, 10.7, and 18.7 GHz. An example of this optimization for DYE-2 in 2021 is presented in Supplementary Fig. 1b.

As a summary, we optimize for each year and each of our 19 sites (Fig. 1) first the number of ice lenses $N_{ice}$ to increase the match between modeled and observed winter $T_B$ at 1.4 GHz, and then the grain size adjustment factor $C_{GS}$ to increase the match between modeled and observed winter $T_B$ at 6.9, 10.7, and 18.7 GHz (Fig. 2). Then for the rest of the year, we extract at time of microwave observations (6 AM) the profiles of snow density, temperature, ice content, and grain size from the GEUS model, process them with these optimal $N_{ice}$ and $C_{GS}$ values and provide the resulting profiles as input to SMRT to get our optimized estimation of $T_B$ (Fig. 2a).

## 3.4. The random forest model ensemble

We use Random Forest (RF) models (Breiman, 2001), an intermediate-complexity machine learning algorithm known for its robustness and ability to handle small, heterogeneous datasets. RF models have been successfully applied to snow depth estimation (Yang et al., 2020) and snow density retrieval (Li et al., 2025). Their relatively low risk of overfitting makes them suitable for our training dataset: the catalogue of corresponding UDLW and $T_B$ produced by our model suite at our 19 sites for the period 2010–2023. To further reduce overfitting and estimate our retrieval's predictive uncertainty, we adopt an ensemble approach (Altman and Krzywinski, 2017): 14 RF models are trained, each excluding one year from the training data. For a given set of input $T_B$, the RF model ensemble provides the mean predicted UDLW and the standard deviation across models, which we interpret as an estimate of uncertainty on unseen data (Fig. 3).

The input variables of each RF model must be carefully selected to help the model extract the relevant information from the training data. To extract the signature of liquid water from the observed $T_B$, many melt detection algorithms use the difference between the measured $T_B$ and
<figure><img src="image_5.png" /><figcaption>Fig. 2. a) Overview of the coupling of the GEUS snow model and SMRT model and use of the tuning parameters N<sub>ice</sub> and C<sub>GS</sub>. b) Conversion of the GEUS model output into an input that can be used by SMRT.</figcaption></figure>

the average T<sub>B</sub> over the preceding winter which serves as a reference for dry conditions (e.g., Zwally and Fiegles, 1994; Torinesi et al., 2003; Leduc-Leballeur et al., 2020; Mousavi et al., 2021). Following a similar approach, we provide to each RF model the V-pol T<sub>B</sub> but also T<sub>B, norm</sub>, the T<sub>B</sub> normalized by the previous winter's T<sub>B</sub>:

$$T_{B,norm}(f) = (T_{B,V-pol}(f) - T_{B,V-pol,winter}(f)) / (273.15 - T_{B,V-pol,winter}(f)) \quad (1)$$

where for a frequency f, T<sub>B, V-pol, winter</sub> is the median of T<sub>B, V-pol</sub> for the first 60 days of that year, and 273.15 K represents the theoretical maximum T<sub>B</sub> of a black body at 0 °C. Each RF model is therefore trained to predict UDLW from 8 inputs: the values of T<sub>B</sub> and T<sub>B, norm</sub> at V-pol for each frequency 1.4, 6.9, 10.7, and 18.7 GHz. When T<sub>B</sub> values are taken from satellite observations, a slight adjustment was needed for the year 2012. AMSR-E failed in October 2011 and AMSR2 was operational only from July 2012 onward, meaning that there is no winter T<sub>B</sub> at 6.9, 10.7,

and 18.7 GHz for that year. To circumvent this data gap, T<sub>B, V-pol, winter</sub> for 2012 is calculated using the last 60 days of 2012.

Our site selection and the higher frequency of short-lived melt episodes imply that lower values of UDLW are over-represented in the training data. This means that the RF models during their training try to fit these ranges of UDLW at the expense of other, less represented UDLW values. To avoid this, we define weights for each training sample, so that all values of UDLW are represented equally in the training set. All samples are sorted based on UDLW values into 0.5 m-thick bins and the weight of each sample is defined as inversely proportional to the number of samples in that bin. The RF models are built with the scikit-learn python package (Pedregosa et al., 2011) and, after testing more than 190 combinations of hyperparameters using 2021 as a hold out, the following design was used: each RF model is composed of 100 trees, with a maximum depth of 15, and a minimum of 20 samples in the final nodes.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Overview of the our retrieval: the random forest (RF) model ensemble, its training data and procedure, and how it is used to estimate the upper depth of liquid water (UDLW) on the ice sheet.</figcaption></figure>

The evaluation of the RF models and their ensemble is done in two different ways. First, each RF model can be evaluated using the UDLW and T<sub>B</sub> produced by the GEUS and SMRT models from their respective hold-out years. This is done by computing the Root Mean Squared Error (RMSE) and the Mean Error (ME):

$$RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (UDLW_{i}^{RF} - UDLW_{i}^{GEUS})^{2}}, \quad (2)$$

$$ME = \frac{1}{n} \sum_{i=1}^{n} (UDLW_{i}^{RF} - UDLW_{i}^{GEUS}), \quad (3)$$

where UDLW<sub>i</sub><sup>GEUS</sup> are the “true” UDLW values in the hold-out year, directly extracted from the GEUS model output and UDLW<sub>i</sub><sup>RF</sup> is the UDLW value retrieved by the RF model from the corresponding T<sub>B</sub> and T<sub>B, norm</sub> values, which SMRT estimated from the same GEUS model output. This gives an overview of the performance of the RF models when given as input the T<sub>B</sub> and T<sub>B, norm</sub> calculated by the same model chain that produced the training data. In a second step, we evaluate the performance of the RF model ensemble when given the observed T<sub>B</sub> and T<sub>B, norm</sub> as input at three training sites where subsurface temperatures were available and at one site which was not used for training (see section 2.4.).

# 4. Results

## 4.1. Evaluation of the RF models on their hold-out year

The performance of the RF models on their respective hold-out sets provides the most representative estimate of the uncertainty expected when applying the models to previously unseen inputs. Table 1 presents the RMSE and ME for all 14 RF models, evaluated on their hold-out sets for both all samples and those with UDLW <5 m. As an example, the scatterplots corresponding to the hold-out year 2019 are presented in Supplementary Figs. 2. For all samples, RMSE ranges from 0.57 to 5.00 m and absolute ME from 0.04 to 1.31 m. These values decrease to 0.54–1.32 m and 0.02–0.37 m, respectively, when considering only samples with UDLW <5 m. This suggests that the retrieval struggles to detect liquid water below a certain depth, likely due to the limited

sensitivity of microwave emission when the water is below a certain depth (see Section 5.2). We also find that the higher RMSE and ME calculated on all samples for some hold-out years, for instance between 2018 and 2023, correspond to the presence in those years of very deep water (UDLW >10 m) which the RF models typically fail to reproduce, causing poor performance statistics. Contrastingly, RMSE and ME values are less variable for samples with UDLW <5 m. We further examine how the prediction error depends on UDLW and on the total liquid water in the simulated profile (Supplementary Fig. 3). The retrieval mean absolute error increases by a factor 8 when UDLW exceeds 5 m compared to shallower conditions, and likewise increases by a factor 5 when the total liquid water is below 0.1 mm compared to profiles with higher water content. This confirms the retrieval's limitations for detecting either very deep water and very low liquid water amount.

## 4.2. Evaluation of the RF model ensemble at four selected sites

FA-13 and FA-15-1 experience ~0.5 m w.e. of summer melt according to the GEUS snow model output, also visible from the long

<table><caption>**Table 1**
Root Mean Squared Error (RMSE) and Mean Error (ME) of the 14 RF models when evaluated on their hold-out set.</caption><thead><tr><th rowspan="2">Hold-out year</th><th colspan="2">All samples</th><th colspan="2">Samples with UDLW &lt;5 m</th></tr><tr><th>RMSE (m)</th><th>ME (m)</th><th>RMSE (m)</th><th>ME (m)</th></tr></thead><tbody><tr><td>2010</td><td>0.57</td><td>0.13</td><td>0.70</td><td>0.29</td></tr><tr><td>2011</td><td>0.84</td><td>0.42</td><td>0.68</td><td>0.31</td></tr><tr><td>2012</td><td>1.38</td><td>0.41</td><td>0.61</td><td>0.29</td></tr><tr><td>2013</td><td>1.36</td><td>0.24</td><td>0.86</td><td>-0.02</td></tr><tr><td>2014</td><td>1.06</td><td>0.13</td><td>0.65</td><td>-0.08</td></tr><tr><td>2015</td><td>1.26</td><td>0.32</td><td>0.92</td><td>0.23</td></tr><tr><td>2016</td><td>1.78</td><td>0.57</td><td>0.54</td><td>0.22</td></tr><tr><td>2017</td><td>0.79</td><td>-0.22</td><td>0.59</td><td>-0.09</td></tr><tr><td>2018</td><td>1.87</td><td>0.04</td><td>0.55</td><td>-0.18</td></tr><tr><td>2019</td><td>3.12</td><td>-1.19</td><td>1.32</td><td>-0.20</td></tr><tr><td>2020</td><td>2.14</td><td>0.12</td><td>1.10</td><td>-0.37</td></tr><tr><td>2021</td><td>4.85</td><td>-0.14</td><td>0.79</td><td>-0.10</td></tr><tr><td>2022</td><td>3.80</td><td>-0.97</td><td>0.61</td><td>-0.09</td></tr><tr><td>2023</td><td>5.00</td><td>-1.31</td><td>0.68</td><td>-0.09</td></tr><tr><td>Median value:</td><td>2.01</td><td>-0.05</td><td>0.68</td><td>-0.09</td></tr></tbody></table>
period of high T<sub>B</sub> (Fig. 4a, c). Thermistor strings installed at both sites show that a large portion of the snow becomes temperate due to melt-water infiltration (Fig. 4b,d colored dots). The binary melt indicators for 18.7, 6.9 and 1.4 GHz are presented as colored dots in Fig. 4a,c and the indicator at 1.4 GHz shows where UDLW can be retrieved. At FA-13, the retrieved UDLW indicates liquid water close to the surface from mid-June to beginning of August, in agreement with the surface wetness status derived at 18.7 GHz (surface wetness at 6.9 GHz was not available for that pixel). From August 15, surface melt ceases and the retrieved UDLW increases from less than 1 m to 4 m in mid-September. These UDLW values agree well with the subsurface temperature measurements, which also show temperatures above -0.2 °C down to 3.5 m depth in mid-September. Past that date, the 1.4 GHz surface wetness indicator (blue dots in Fig. 4a) does not detect liquid water and we can see that the retrieved values, even though discarded, fail to follow the further increase of UDLW to greater depth (gray dots in Fig. 4b).

At FA-15-1, our retrieval indicates liquid water close to the surface, with UDLW<0.5 m from the beginning of the melting season in mid-June to the beginning of August. Then, surface melt stopped for three weeks, as indicated by the surface wetness status at 18.7 and 6.9 GHz and lower T<sub>B</sub>. During that period, both our retrieved UDLW and the thermistor string indicate refreezing close to the surface and the presence of liquid water from 1.5 to 2 m depth. From the 29th of August to the 2nd of September, another brief melt event occurred, as indicated by the surface wetness status at 18.7 and 6.9 GHz and thermistor data. Potentially due to the smoothing of the T<sub>B</sub>, the retrieved UDLW fails to capture the short-lived surface meltwater and only shows a slightly shallower wet layer, before the surface melt stops and deeper refreezing water is again detected by our retrieval. Interestingly, this short melt event appears to last longer in the thermistor data compared to the binary surface wetness indicator at 18.7 and 6.9 GHz. This could reflect the heterogeneity of surface melt in the area and explain the deviation between the refreezing front as shown by the thermistor data and our retrieved UDLW.

We now evaluate our retrieval at a site receiving significantly less melt than FA-13 and FA-15-1: CP1 for the summer of 2021 (Fig. 5). During that summer, several short-lived melt events occurred in June and July, as visible from the T<sub>B</sub> peaks at most frequencies (Fig. 5a). In mid-August, a more pervasive melt event occurred, generating enough meltwater for some to remain liquid for several weeks, as indicated by the gradual T<sub>B</sub> decrease at all frequencies (Fig. 5a). The thermistors indicate significant melt, with isothermal conditions (interpreted as infiltration) down to 3.5 m in mid-August. After the melt peak, the thermistor string shows a reduction in the isothermal portion of the firn, with cooling/refreezing at its upper and lower boundaries. But unlike other instruments, the thermistor at a depth of 2.5 m still indicated 0 °C at the beginning of September, indicating remnant subsurface liquid water.

Our retrieved UDLW successfully captured near-surface meltwater during the first melt event of the season. During the second event, however, the retrieval was challenged, as many estimates had too high a standard deviation and were discarded (gray dots). For the third and slightly longer event in mid-August, the retrieval indicated shallow water (UDLW <1.5 m) instead of surface water, as would normally be expected during surface melt. Despite the binary surface wetness indicator at 18.7 GHz showing continuous surface water between 14 and 22 August, the 18.7 GHz T<sub>B</sub> (Fig. 5a) showed a more complex pattern: T<sub>B</sub> > 265 K on 14–15 August, slightly lower T<sub>B</sub> on 16–18 August, a brief return above 265 K on the 19th, and then a sharp drop below 248 K from the 20th onward. This illustrates how binary thresholding can cause information loss and inaccurate timing compared to a detailed analysis of the T<sub>B</sub> time series. During that event, the retrieved UDLW closely mirrored the T<sub>B</sub> evolution and provided a more nuanced picture than the binary indicator. After the melt event, the UDLW increased first to 1.6 m and then to 3 m, indicating gradual refreezing of meltwater in the upper snow layer and exposing the underlying wet layer. The dry surface

conditions indicated by the 18.7 GHz T<sub>B</sub> on 20–26 August contrast with the thermistor data which showed isothermal, potentially wet, conditions close to the surface. This could stem from a local effect around the thermistor string, which is not representative of the satellite footprint. Yet, the thermistor data, along with the surface wetness indicators at 6.9 and 1.4 GHz, confirm the presence of subsurface water and show that our algorithm can provide meaningful estimates of UDLW.

After validating our retrieval against independent measurements at three training sites, we examine its performance during the 2021 melt season at the Qassimiut Aquifer location (Fig. 1), which was not among the sites used for the training of the RF model ensemble. No in situ temperature or water-content data exist there, so we compare the retrieved UDLW with liquid-water profiles simulated by the GEUS snow model (Fig. 6). During the active melt period (20 June–1 September), the retrieval places the UDLW within 0.6 m of the surface, consistent with the model, which also shows near-surface liquid water. After melt ceases, both datasets show a recession of the wetted layer similar to that observed at FA-13 and FA-15-1 (Fig. 4). The retrieval, however, indicates a faster increase in UDLW than the model. This difference may reflect model shortcomings, such as insufficient upward heat conduction or too much retained meltwater, both of which would slow refreezing. The retrieval becomes unreliable below ~4 m and therefore cannot detect the upper boundary of the aquifer simulated at 7 m depth. Given uncertainties in the model's infiltration scheme and the lack of lateral flow, the simulated water content remains highly uncertain. In this setting the retrieval provides new insight: meltwater infiltrates to at least 4 m depth, and the rapid post-melt recession of the wetted layer suggests either efficient heat loss to the surface or effective drainage that leaves little water to refreeze in autumn. This highlights the retrieval's realism beyond the training sites and its value for evaluating reanalysis-driven snow model simulations.

### 4.3. Evolution of UDLW on the Greenland ice sheet

The RF model ensemble is then applied between 2010 and 2023 to all pixels in the accumulation area of the Greenland Ice Sheet. Fig. 7 shows, for each year, the number of days with estimated UDLW deeper than 1 m, reflecting both the melt intensity and the infiltration depth for that year. As expected, high melt years such as 2010, 2012, 2019, and 2021 exhibit a larger number of days with deep water. Note that in 2012, measurements from AMSR and consequently our retrieval are only available from July onwards. The southeastern region consistently shows long periods where water below 1 m is detected, consistent with observations in that region of meltwater infiltration down to 12 m to sustain firn aquifers (Miller et al., 2020). Deep percolation is also visible in the western percolation area which spans ~1800–2500 m a.s.l. in the southern half of the ice sheet and between 1300 and 1700 m a.s.l. in the northwest. This deep percolation area is bounded at higher elevations by limited meltwater supply. We also note a fringe at the lower boundary of the accumulation area where SMOS, at 1.4 GHz, does not detect water in any year. This is due to a change in how the passive microwave signal responds to surface melt in these rougher, more ice-dominated areas of the Greenland Ice Sheet (Miller et al., 2022), which makes the wetness indicator unable to detect melt.

The retrieved UDLW dataset can also be used to study the seasonal evolution of percolation and refreezing. Here, we use three dates from 2021, the 16th, 20th and 25th of August, as an illustration of our retrieval capacity (Fig. 8). On the 16th of August, liquid water was detected in most of the percolation zone in western Greenland and on all of the southern dome (Fig. 8a). The depth of liquid water, as estimated by our RF model ensemble, ranges between 0 m where the melt is more intense to 1 m in the uppermost melting regions. Four days later, on the 20th of August 2021, our retrieval indicates that only the lower part of the percolation zone still has water within 1 m of the surface (in yellow in Fig. 8b). The upper percolation area of the ice sheet shows liquid water at about 2 m depth, indicated by dark orange in Fig. 8b. The
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Evaluation of the UDLW retrieval at FA-13 in 2013 (a, b) and FA-15-1 in 2015 (c, d). Panels a and c: Brightness temperatures ($T_B$) at 1.4, 6.9, 10.7, and 18.7 GHz before smoothing (light colored lines) and after smoothing (deep colored lines), with binary melt indicators on top. Panels b and d: Upper depth of the liquid water (UDLW) estimated by our RF model ensemble (dots for the mean, gray shading for the intra-ensemble standard deviation), when 1.4 GHz indicates a wet surface (black dots) or dry surface (gray dots), along with the measured subsurface temperature (colored dots) in the background, red indicating temperatures above –0.2 °C, indicative of infiltrating or refreezing meltwater. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Evaluation of the UDLW retrieval at CPI in 2021, a) Brightness temperatures ($T_B$) at 1.4, 6.9, 10.7, and 18.7 GHz before smoothing (light colored lines) and after smoothing (deep colored lines). b) Upper depth of the liquid water (UDLW) estimated by our RF model ensemble (dots for the mean, gray shading for the intra-ensemble standard deviation), when 1.4 GHz indicates a wet surface (black dots) or dry surface (gray dots), along with the measured subsurface temperature (colored dots) in the background, red indicating temperatures above -0.2 °C, indicative of infiltrating or refreezing meltwater. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Evaluation of the UDLW retrieval at the Qassimiut Aquifer site in 2021, a) Brightness temperatures ($T_B$) at 1.4, 6.9, 10.7, and 18.7 GHz before smoothing (light colored lines) and after smoothing (deep colored lines). b) Upper depth of the liquid water (UDLW) estimated by our RF model ensemble (dots for the mean, gray shading for the intra-ensemble standard deviation), when 1.4 GHz indicates a wet surface (black dots) or dry surface (gray dots), along with the liquid water content simulated by the CARRA-forced GEUS snow model.</figcaption></figure>
summit of the southern dome has large retrieval uncertainties, shown in dark gray in Fig. 8b and in green in Fig. 8e, but is still classified as wet by the SMOS binary indicator. For regions with UDLW greater than 1 m, our retrieval algorithm confirms that surface melt has likely stopped and that only refreezing subsurface water is left. Five days later, on the 25th of August (Fig. 8c), the area where liquid water is detected has shrunk to a lower elevation, and the summit of the southern dome of the ice sheet is now categorized as dry. The retrieved UDLW in the remaining areas with liquid water is deeper than 2 m for most pixels. In the upper percolation area, there is a band of pixels, shown in dark gray in Fig. 8c and in yellow in Fig. 8f, that is classified as wet by the SMOS binary indicator and for which our retrieval has too large uncertainties (>1 m) to provide a UDLW value. We conclude that on that day, most of the pixels marked as wet by the SMOS binary indicator, are, in fact, infiltrated subsurface liquid water and not surface melt.

# 5. Discussion

We showed that by using a catalogue of simulations from coupled snow and radiative transfer models, it is possible to train an algorithm capable of estimating UDLW from observed multi-frequency T<sub>B</sub>. Our work represents an improvement compared to previous retrievals as it gives a quantitative estimate of the UDLW instead of categorical descriptions of subsurface water as provided by Colliander et al. (2023) and Leduc-Leballeur et al. (2025). Additionally, the use of a model suite to build a diverse training set makes our algorithm applicable to the entire ice sheet, representing a significant improvement over Shang et al. (2025), whose algorithm was trained specifically to reproduce seven years of airborne radar data from the aquifer region in Greenland.

In the remainder of this section, we first discuss the relative importance of the input variables for the RF model ensemble and how optimizing the coupling between SMRT and the GEUS snow model improves the training catalogue. We then assess the depth range to which our retrieval, and the passive microwave observations it is based on, is sensitive. Lastly, we outline future work that could refine or expand the retrieval approach.

## 5.1. Feature importance in the RF models and relevance of the optimization of the coupling between SMRT and the GEUS snow model

Random Forest models offer a straightforward way to assess the importance of each input variable. Importance is defined as the total reduction in mean squared error (MSE) that a variable contributes across all the splits where it is used during training, summed over all trees in the ensemble. The variable that leads to the largest overall reduction in MSE is considered the most important. The most important variables in our RF model ensemble are the absolute and normalized 1.4 GHz T<sub>B</sub> which are responsible, on average across RF models, for 65 % of the RF model ensemble performance. The absolute T<sub>B</sub> at 10.7 and 6.9 GHz explain 18 % of the performance together, and the remaining variables bring minor contributions between 2 and 6 %. This ranking of the input variables reflects our focus on the detection of subsurface water and that the 1.4 GHz T<sub>B</sub> is the input variable that is most sensitive to deep water (e.g. Picard et al., 2022a). This is also another reminder that the predicted UDLW is tightly linked to the absolute and normalized values of T<sub>B</sub> that are given as input to the RF model ensemble, and therefore that our retrieval will perform best if its training T<sub>B</sub> are similar to the observed T<sub>B</sub> on which it will be applied.

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Number of days with estimated subsurface liquid water depth (UDLW) deeper than 1 m each year from 2010 to 2023.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Maps of the depth of subsurface liquid water (UDLW) estimated by the RF model ensemble (a-c) and associated uncertainty (d-f) for three days in August 2021.</figcaption></figure>

We improve the similarity between observed and simulated T<sub>B</sub>
through the optimization of our coupling between SMRT and the GEUS
snow model. The gain of the optimization is illustrated in Fig. 9 where T<sub>B</sub>
simulated without and with optimization of N<sub>ice</sub> and C<sub>GS</sub> are compared
to satellite observations at four sites: KAN_U, DYE-2, CP1 and Saddle
(Fig. 1 for locations). All histograms show a bi-modal distribution, with
the first mode corresponding to low T<sub>B</sub> occurring in winter and the
second mode corresponding to high T<sub>B</sub> during the melting season.
Without optimization, when SMRT uses the imperfect grain size profile
from the GEUS snow model output and assumes a single ice layer per
depth interval, the simulated winter T<sub>B</sub> values show a systematic offset
compared to the observations (low-value modes in orange histograms
compared to the black ones in Fig. 9). We found that RF models trained
on these simulated T<sub>B</sub> were producing nonsensical UDLW values when
given observed T<sub>B</sub> as input. Contrasting with the non-optimized T<sub>B</sub>, the
T<sub>B</sub> values produced when optimizing N<sub>ice</sub> and C<sub>GS</sub> are visibly closer to
the observed T<sub>B</sub> (low-value modes in green histograms compared to the
black ones in Fig. 9). RF models trained on these optimized simulated T<sub>B</sub>
will, when used on observed T<sub>B</sub>, produce UDLW values that are closer to
the realistic UDLW values contained in the training catalogue.

5.2. Sensitivity of the T<sub>B</sub> observations to the deep water

Our retrieval shows its best performance when the upper bound of
the wetted layer is within 5 m of the surface, both when evaluated
against model-produced data (Table 1) and observed T<sub>B</sub> (Fig. 4). This
prompts the question of the depth range for which multi-frequency
passive microwave data can be used to detect subsurface liquid water.
The maximum depth to which our retrieval is possible can be explored,
as in Colliander et al. (2022), by comparing observed T<sub>B</sub> with the depth

of the dry-to-wet transition derived from thermistor strings (assuming
that temperatures higher than -0.2 °C indicate liquid water). This
analysis is possible at FA-13 and FA-15-1 (Fig. 4), where deep and
vertically dense thermistor strings allow us to track the progression of
the refreezing front after the melt season ends (Fig. 4b, d), from the
surface to ~12 m, where the perennial firn aquifer lies (e.g., Miller et al.,
2020). Observed T<sub>B</sub> decrease at different rates depending on frequency
as the dry-to-wet transition deepens (Fig. 4): T<sub>B</sub> at 18.7 GHz drops
sharply once the transition recedes below 1 m, while T<sub>B</sub> at 10.7, 6.9, and
1.4 GHz decreases more gradually, reflecting the deepening wet front
and declining liquid water content. Notably, the rates of T<sub>B</sub> decrease at
these three frequencies are rather similar (Fig. 4a, c) although these
three frequencies are expected, in theory, to have different sensitivity to
the depth and amount of buried liquid water (Picard et al., 2022a). By
December 1st at FA-13 and February 1st at FA-15-1, T<sub>B</sub> at 10.7, 6.9, and
1.4 GHz stabilizes around ~230 K. In contrast, temperature measure-
ments show the dry-to-wet front continuing to descend past 7 m. The
further retreat of the wet front beyond ~7 m has no apparent impact on
T<sub>B</sub> at these frequencies. These observations suggest that, at FA-13 and
FA-15-1, the maximum detection depth for subsurface meltwater using
T<sub>B</sub> at 10.7, 6.9, and 1.4 GHz is approximately 7 m; meltwater deeper
than this leaves no clear radiometric signature at these frequencies.

5.3. *Other limitations and future work*

Despite our efforts to maximize the similarity between simulated and
observed T<sub>B</sub>, Fig. 9 shows that some mismatch remains. Discrepancies
between model and observations represent the main limiting factor
when training a retrieval on model-produced data. Therefore future
work should aim to reduce this discrepancy by investigating additional
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Histograms of brightness temperatures ($T_B$) observed by SMOS and AMSR, simulated by the model suite without optimization of $N_{\text{ice}}$ and $C_{\text{GS}}$, and simulated with optimization of $N_{\text{ice}}$ and $C_{\text{GS}}$ to minimize the mismatch with winter $T_B$ observations.</figcaption></figure>

physical and electromagnetic snowpack characteristics that may be misrepresented in SMRT and the GEUS snow model. For example, ice layers are currently simulated in SMRT as horizontally continuous, whereas several studies have shown that ice disks and pipes can alter the scattering properties of the snowpack (e.g. Zwally, 1977; Jezek et al., 1994). Accounting for such features in the emission modeling could help improve the simulations and provide a better understanding of the co-variation between multi-frequency $T_B$ and subsurface water. Yet, no comprehensive information is currently available on the number, depth distribution, and geometry of these ice inclusions. Another cause of the remaining mismatch between simulated and observed $T_B$ stems from the spatial scales involved. SMRT simulates $T_B$ based on a single evolving uni-dimensional profile, while observed $T_B$ integrates signals over a large footprint (32 km at 1.4 GHz) including significant altitudinal variations. As a result, changes in the modeled profiles, such as the refreezing of all liquid water, can cause abrupt shifts in simulated $T_B$, whereas the observed $T_B$ from a drying pixel will transition more gradually due to sub-pixel heterogeneity. Another impact of the satellites' resolution is the scale gap between the in situ measurements and the passive microwave observations which means that some ground observations are not fully representative of the larger pixel they are located in. This scale gap could also potentially explain the constant

offset between our retrieved UDLW and the observed refreezing front seen at FA-15-1 (Fig. 4d), if the site is in a more melt-prone area than its surroundings. Bridging this resolution gap would require incorporating sub-pixel meltwater dynamics—potentially tied to surface topography—into the radiative transfer model, but current knowledge on their frequency-dependent impact remains limited.

Current evaluations of remote sensing retrievals rely on surface energy balance and snow models forced by regional climate models or reanalysis data (Colosio et al., 2021; Houtz et al., 2021; Mousavi et al., 2021; Hossan et al., 2025). These models are limited by uncertainties in meteorological forcing and the conversion of weather data into melt estimates (e.g., Covi et al., 2023), and their predictions of liquid water content remain poorly constrained due to infiltration scheme variability (Vandecrux et al., 2020b). Measurements of the liquid water profiles (e.g. Samimi et al., 2020) and percolation depth (e.g. Heilig et al., 2018) on the Greenland Ice Sheet would be extremely useful for the evaluation and calibration of both remote sensing retrievals and snow models, but they remain very scarce to this day. Here, we used subsurface temperature observations from thermistor strings—widely deployed at Greenland weather stations—to infer percolation depth and verify that the GEUS snow model produced realistic infiltration patterns at our sites. These observations, although very informative, can also include
measurement errors such as when the meltwater infiltrates along the thermistor string (e.g. Humphrey et al., 2012) or occurrences when thermistors remained below freezing point even though a nearby weather station and passive microwave binary indicators indicate extended periods of melt (e.g. in the 2019 data at CP1 presented in Saito et al., 2024). These cases that needed to be discarded, in combination with thermistors that are too far from the surface, explain the small number of sites and years where evaluation of our retrieval was possible. Improving and expanding the use of thermistors, or reinterpreting existing datasets (e.g., Moon et al., 2025), could improve model evaluation.

We used a Random Forest algorithm as a fast and efficient way to “invert” SMRT, estimating subsurface water depth from four vertically polarized T<sub>B</sub> and their normalized counterparts (eight predictors). In the future, full inversion of SMRT may be explored but it would require resolving the liquid water profile across multiple layers while also accounting for snow grain size, ice layer characteristics, and ice inclinations—parameters that must be jointly tuned to reproduce observed multi-frequency emissions in percolation zones. This complexity could be addressed through Bayesian inversion using firn model outputs as priors (e.g., Singh et al., 2024). Another potential development would be to use H-pol T<sub>B</sub> or microwave measurements at higher frequencies (e.g., 36.5 and 89.0 GHz). Yet, H-pol T<sub>B</sub> are more sensitive to snow layering, and higher frequencies are strongly influenced by snow microstructure and atmospheric conditions (Picard et al., 2022a; Leduc-Leballeur et al., 2025), meaning that their inclusion in a retrieval technique would require careful treatment of these additional sources of variability, with an unknown cost-to-benefit ratio.

## 6. Conclusion

Meltwater infiltration and refreezing beneath the snow surface have a direct influence on runoff from the Greenland Ice Sheet accumulation area and is notoriously difficult to observe at large scale. Here, we combine the GEUS snow model with SMRT to generate a catalogue of simulated snow profiles containing liquid water, from which we derive the upper depth of liquid water (UDLW) and simulate T<sub>B</sub> at 1.4, 6.9, 10.7, and 18.7 GHz. Using this catalogue, we train an ensemble of cross-validated random forest models to estimate the UDLW from V-pol T<sub>B</sub> at four frequencies. Applying our retrieval on observed T<sub>B</sub> from SMOS, AMSR-E, and AMSR2 provides daily UDLW maps from 2010 to 2023 and thereby the first quantitative retrieval of UDLW across the Greenland Ice Sheet accumulation zone. The retrieval agrees well with subsurface temperature records at three selected sites and years. Although restricted to the accumulation area and showing best performance when the water is within five meters of the surface, the daily UDLW maps capture seasonal infiltration and refreezing patterns, reveal new aspects of meltwater storage, and offer a basis for evaluating how meltwater storage influences the surface mass balance.

## Data and code availability

The daily maps of UDLW for 2010–2023 are available at 10.22008/FK2/69CRRT (Vandecrux, 2025a). The simulated T<sub>B</sub> at 19 sites produced by the coupled and optimized SMRT and GEUS snow models are available at 10.22008/FK2/KKIOZZ (Vandecrux, 2025b). The SMOS T<sub>B</sub> and snow wetness grids at enhanced resolution are available at 10.57932/f72f9515-3699-4fae-92fc-a350075d042f (Zeiger and Picard, 2024a) and 10.57932/1970fb7c-cdb4-4ddf-9891-1e6836a46f25 (Zeiger and Picard, 2024b). The original AMSR2 and AMSR-E data are available at https://gportal.jaxa.jp/gpr/ and https://doi.org/10.5067/AMSR-E/AE_L2A.003 (Ashcroft and Wentz, 2013), but we here used the automatically updated, re-gridded combination of these two datasets, available as a zarr file system: https://snow.univ-grenoble-alpes.fr/ope ndata/PMW-Greenland-12km.zarr (Picard, 2024), and downloadable through linux wget or python xarray.open_zarr. The SUMup collaborative

database is available at 10.18739/A2M61BR5M (Vandecrux et al., 2024a). The CARRA reanalysis is available at 10.24381/cds.713858f6 (Schyberg et al., 2020).

The script for training of the RF models and for producing the UDLW maps are available at https://doi.org/10.5281/zenodo.15791029 (Vandecrux and Picard, 2025). The GEUS snow model code is hosted at https://github.com/GEUS-Glaciology-and-Climate/GEUS-SEB-firn-model and the SMRT model at https://github.com/smrt-model/smrt. This research is based on open data, and its results are likewise published in an open data repository.

## CRediT authorship contribution statement

**Baptiste Vandecrux:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. **Ghislain Picard:** Writing – review & editing, Software, Methodology, Data curation, Conceptualization. **Pierre Zeiger:** Writing – review & editing, Data curation. **Marion Leduc-Leballeur:** Writing – review & editing, Data curation. **Andreas Colliander:** Writing – review & editing. **Alamgir Hossan:** Writing – review & editing. **Andreas Ahlstrøm:** Resources, Funding acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This study was supported by the European Space Agency (ESA) Climate Change Initiative (CCI) under the Water Under Snow Cover CCI research fellowship. Support was also received from the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) and Greenland Climate Network (GC-Net) at the Geological Survey of Denmark and Greenland. A contribution to this work was made at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the USA National Aerospace and Space Administration. We thank Anna Maria Trofaier and Clément Albergel, our technical officers at ESA, and Jason Box for their support. We are grateful to our editor, Menghua Wang, and three anonymous reviewers for their constructive feedback.

## Appendix A. Supplementary data

Supplementary figures to this article can be found online at https://doi.org/10.1016/j.rse.2025.115197.

## References

Altman, N., Krzywinski, M., 2017. Ensemble methods: bagging and random forests. Nat. Methods 14, 933–934. https://doi.org/10.1038/nmeth.4438.

[dataset] Ashcroft, P., Wentz, F.J., 2013. AMSR-E/aqua L2A global swath spatially-resampled brightness temperatures (AE_L2A, version 3). NASA National Snow and Ice Data Center Distributed Active Archive Center. https://doi.org/10.5067/A MSR-E/AE_L2A.003.

Box, J.E., Hubbard, A., Bahr, D.B., Aschwanden, A., Cuzzone, J., Grinsted, A., Hofer, S., Solgaard, A.M., van As, D., van den Broeke, M.R., Wilton, D., 2022. Greenland ice sheet climate disequilibrium and committed sea-level rise. Nat. Clim. Chang. 12, 808–813. https://doi.org/10.1038/s41558-022-01441-2.

Breiman, L., 2001. Random forests. Mach. Learn. 45, 5–32. https://doi.org/10.1023/A:1010933404324.

Brucker, L., Picard, G., Arnaud, L., Barnola, J., Schneebeli, M., Brunjail, H., et al., 2011. Modeling time series of microwave brightness temperature at dome C, Antarctica, using vertically resolved snow temperature and microstructure measurements. J. Glaciol. 57, 171–182. https://doi.org/10.3189/002214311795306736.

Charalampidis, C., Van As, D., Colgan, W.T., Fausto, R.S., Macferrin, M., Machguth, H., 2016. Thermal tracing of retained meltwater in the lower accumulation area of the southwestern Greenland ice sheet. Ann. Glaciol. 57, 1–10. https://doi.org/10.1017/aog.2016.2.
Chen, N., Li, W., Fan, Y., Zhou, Y., Aoki, T., Tanikawa, T., Niwano, M., Hori, M., Shimada, R., Matoba, S., Stamnes, K., 2025. Snow parameter retrieval (SPR) algorithm for the GCOM-S/C-GSLI sensor: validation over the Greenland ice sheet. Front. Environ. Sci. 13, 1541041. https://doi.org/10.3389/fenvs.2025.1541041.

Colliander, A., Mousavi, M., Marshall, S., Samimi, S., Kimball, J.S., Miller, J.Z., Johnson, J., Burgin, M., 2022. Ice sheet surface and subsurface melt water discrimination using multi-frequency microwave radiometry. Geophys. Res. Lett. 49. https://doi.org/10.1029/2021GL096599 e2021GL096599.

Colliander, A., Mousavi, M., Kimball, J.S., Miller, J.Z., Burgin, M., 2023. Spatial and temporal differences in surface and subsurface meltwater distribution over Greenland ice sheet using multi-frequency passive microwave observations. Remote Sens. Environ. 295, 113705. https://doi.org/10.1016/j.rse.2023.113705.

Colosio, P., Tedesco, M., Fettweis, X., Ranzi, R., 2021. Surface melting over the Greenland ice sheet derived from enhanced resolution passive microwave brightness temperatures (1979–2019). Cryosphere 15, 2623–2646. https://doi.org/10.5194/tc-15-2623-2021.

Covi, F., Hock, R., Reijmer, C.H., 2023. Challenges in modeling the energy balance and melt in the percolation zone of the Greenland ice sheet. J. Glaciol. 69, 164–178. https://doi.org/10.1017/jog.2022.67.

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.O., Shields, C.L., Solgaard, A.M., Box, J.E., 2021. Programme for monitoring of the Greenland ice sheet (PROMICE) automatic weather station data. Earth Syst. Sci. Data 13, 3819–3845. https://doi.org/10.5194/essd-13-3819-2021.

Heilig, A., Eisen, O., MacFerrin, M., Tedesco, M., Fettweis, X., 2018. Seasonal monitoring of melt and accumulation within the deep percolation zone of the Greenland ice sheet and comparison with simulations of regional climate modeling. Cryosphere 12, 1851–1866. https://doi.org/10.5194/tc-12-1851-2018.

Horlings, A.N., Christianson, K., Miège, C., 2022. Expansion of firn aquifers in Southeast Greenland. Journal of Geophysical Research: Earth Surface 127. https://doi.org/10.1029/2022JF006753 e2022JF006753.

Hossan, A., Colliander, A., Vandecrux, B., Schlegel, N.-J., Harper, J., Marshall, S., Miller, J.Z., 2025. Retrieval and validation of total seasonal liquid water amounts in the percolation zone of Greenland ice sheet using L-band radiometry. Cryosphere 19, 4237–4258. https://doi.org/10.5194/tc-19-4237-2025.

Houtz, D., Mätzler, C., Naderpour, R., Schwank, M., Steffen, K., 2021. Quantifying surface melt and liquid water on the Greenland ice sheet using L-band radiometry. Remote Sens. Environ. 256, 112341. https://doi.org/10.1016/j.rse.2021.112341.

Humphrey, N.F., Harper, J.T., Pfeffer, W.T., 2012. Thermal tracking of meltwater retention in Greenland's accumulation area. Journal of Geophysical Research: Earth Surface 117, F01010. https://doi.org/10.1029/2011JF002083.

Jezek, K.C., Gogineni, P., Shanableh, M., 1994. Radar measurements of melt zones on the Greenland ice sheet. Geophys. Res. Lett. 21, 33–36. https://doi.org/10.1029/93GL03377.

Kerr, Y.H., Waldteufel, P., Wigneron, J.-P., Martinuzzi, J., Font, J., Berger, M., 2001. Soil moisture retrieval from space: the soil moisture and ocean salinity (SMOS) mission. IEEE Trans. Geosci. Remote Sens. 39, 1729–1735. https://doi.org/10.1109/36.942551.

Køltzow, M., Schyberg, H., Støylen, E., Yang, X., 2022. Value of the Copernicus Arctic regional reanalysis (CARRA) in representing near-surface temperature and wind speed in the north-east European Arctic. Polar Res. 41, 8002. https://doi.org/10.33265/polar.v41.8002.

Kontu, A., Pulliainen, J., 2008. Simulation of spaceborne microwave radiometer measurements of snow cover using in situ data and brightness temperature modeling. IEEE Trans. Geosci. Remote Sens. 48, 1031–1044. https://doi.org/10.1109/IGARSS.2008.4780021.

Langen, P.L., Fausto, R.S., Vandecrux, B., Mottram, R.H., Box, J.E., 2017. 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. https://doi.org/10.3389/feart.2016.01010.

Langlois, A., Royer, A., Derksen, C., Montpetit, B., Dupont, F., Goïta, K., 2012. Coupling the snow thermodynamic model SNOWPACK with the microwave emission model of layered snowpacks for subarctic and arctic snow water equivalent retrievals. Water Resour. Res. 48, W12524. https://doi.org/10.1029/2012WR012133.

Leduc-Leballeur, M., Picard, G., Macelloni, G., Arnaud, L., Brogioni, M., Mialon, A., Kerr, Y.H., 2017. Influence of snow surface properties on L-band brightness temperature at dome C. Antarctica. Remote Sens. Environ. 199, 427–436. https://doi.org/10.1016/j.rse.2017.07.026.

Leduc-Leballeur, M., Picard, G., Macelloni, G., Mialon, A., Kerr, Y.H., 2020. Melt in Antarctica derived from soil moisture and ocean salinity (SMOS) observations at L band. Cryosphere 14, 539–548. https://doi.org/10.5194/tc-14-539-2020.

Leduc-Leballeur, M., Picard, G., Zeiger, P., Macelloni, G., 2025. Empirical classification of dry-wet snow status in Antarctica using multi-frequency passive microwave observations. EGUSphere [preprint]. https://doi.org/10.5194/egusphere-2025-732.

Li, W., Veldhuijsen, S.B.M., Lhermitte, S., 2025. Machine learning of Antarctic firn density by combining radiometer and scatterometer remote-sensing data. Cryosphere 19, 37–61. https://doi.org/10.5194/tc-19-37-2025.

Löwe, H., Picard, G., 2015. Microwave scattering coefficient of snow in MEMLS and DMRT-ML revisited: the relevance of sticky hard spheres and tomography-based estimates of stickiness. Cryosphere 9, 2101–2117. https://doi.org/10.5194/tc-9-2101-2015.

MacFerrin, M., Machguth, H., Van As, D., Charalampidis, C., Stevens, C.M., Heilig, A., Vandecrux, B., Langen, P.L., Mottram, R., Fettweis, X., Van den Broeke, M.R., 2019. Rapid expansion of Greenland's low-permeability ice slabs. Nature 573, 403–407. https://doi.org/10.1038/s41586-019-1550-3.

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

Mätzler, C., 1998. Improved born approximation for scattering of radiation in a granular medium. J. Appl. Phys. 83, 6111–6117. https://doi.org/10.1063/1.367496.

Meyer, H., Pebesma, E., 2021. Predicting into unknown space: Estimating the area of applicability of spatial prediction models. Methods Ecol. Evol. 12, 1620–1633. https://doi.org/10.1111/2041-210X.13650.

Miège, C., Forster, R.R., Brucker, L., Koenig, L.S., Solomon, D.K., Paden, J.D., Box, J.E., Burgess, E.W., Miller, J.Z., McNerney, L., Brautigam, N., Fausto, R.S., Gogineni, S., 2016. Spatial extent and temporal variability of Greenland firn aquifers detected by ground and airborne radars. Journal of Geophysical Research: Earth Surface 121, 2381–2398. https://doi.org/10.1002/2016JF003869.

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., McConnell, J.R., 2020. Hydrology of a perennial firn aquifer in Southeast Greenland: an overview driven by field data. Water Resour. Res. 56. https://doi.org/10.1029/2019WR026348 e2019WR026348.

Miller, J.Z., Culberg, R., Long, D.G., Shuman, C.A., Schroeder, D.M., Brodzik, M.J., 2022. An empirical algorithm to map perennial firn aquifers and ice slabs within the Greenland ice sheet using satellite L-band microwave radiometry. Cryosphere 16, 103–125. https://doi.org/10.5194/tc-16-103-2022.

Moon, T., Harper, J., Colliander, A., Hossan, A., Humphrey, N., 2025. L-band radiometric measurement of liquid water in Greenland's firn: Comparative analysis with in situ measurements and modeling. Ann. Glaciol. 66, e15. https://doi.org/10.1017/aog.2025.10012.

Mote, T.L., Anderson, M.R., Kuivinen, K.C., Rowe, C.M., 1993. Passive microwave-derived spatial and temporal variations of summer melt on the Greenland ice sheet. Ann. Glaciol. 17, 233–238. https://doi.org/10.1389/S0260305500012891.

Mousavi, M., Colliander, A., Miller, J.Z., Entekhabi, D., Johnson, J.T., Shuman, C.A., Kimball, J.S., 2021. Evaluation of surface melt on the Greenland ice sheet using SMAP L-band microwave radiometry. IEEE J. Sel. Top. Appl. Earth Res. Remote Sens. 14, 11439–11449. https://doi.org/10.1109/JSTARS.2021.3124229.

Osotaka, I.N., Shepherd, A., Ivins, E.R., Schlegel, N.-J., Amory, C., van den Broeke, M.R., Horwath, M., Joughin, I., King, M.D., Krinner, G., Nowicki, S., Payne, A.J., Rignot, E., Scambos, T., Simon, K.M., Smith, B.E., Sørensen, L.S., Velicogna, I., Whitehouse, P.L., A., Agosta, C., Ahlström, A.P., Blazquez, L.S., Colgan, W., Engdahl, M.E., Fettweis, X., Forsberg, R., Gallée, H., Gardner, A., Gilbert, L., Gournelen, N., Groh, A., Gunter, B.C., Harig, C., Helm, V., Khan, S.A., Kittel, C., Konrad, H., Langen, P.L., Lecavalier, B.S., Liang, C.-C., Loomis, B.D., McMillan, M., Melini, D., Mernild, S.H., Mottram, R., Mouginot, J., Nilsson, J., Noël, B., Pattle, M. E., Peltier, W.R., Pie, N., Rocha, M., Sasgen, I., Save, H.V., Seo, K.-W., Scheuchl, B., Schrama, E.J.O., Schröder, L., Simonsen, S.B., Slater, T., Spada, G., Sutterley, T.C., Vishwakarma, D.B., van Wessem, J.M., Wiese, D., van der Wal, W., Wouters, B., 2023. Mass balance of the Greenland and Antarctic ice sheets from 1992 to 2020. Earth Syst. Sci. Data 15, 1597–1616. https://doi.org/10.5194/essd-15-1597-2023.

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, A., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., 2011. Scikit-learn: machine learning in Python. J. Mach. Learn. Res. 12, 2825–2830. https://doi.org/10.48550/arXiv.1201.0490.

Picard, G., Sandells, M., Löwe, H., 2018. SMRT: an active-passive microwave radiative transfer model for snow with multiple microstructure and scattering formulations (v1.0). Geosci. Model Dev. 11, 2763–2788. https://doi.org/10.5194/gmd-11-2763-2018.

Picard, G., Leduc-Leballeur, M., Banwell, A.F., Brucker, L., Macelloni, G., 2022a. The sensitivity of satellite microwave observations to liquid water in the Antarctic snowpack. Cryosphere 16, 5061–5083. https://doi.org/10.5194/tc-16-5061-2022.

Picard, G., Löwe, H., Domine, F., Arnaud, L., Larue, F., Favier, V., et al., 2022b. The microwave snow grain size: a new concept to predict satellite observations over snow-covered regions. AGU Adv. 3. https://doi.org/10.1029/2021AV000630 e2021AV000630.

Polder, D., Van Santeen, J.H., 1946. The effective permeability of mixtures of solids. Physica 12, 257–271. https://doi.org/10.1016/S0031-8914(46)80066-1.

Saito, J., Harper, J., Humphrey, N., 2024. Uptake and transfer of heat within the firn layer of the Greenland ice sheet's percolation zone. J. Geophys. Res.: Earth Surf. 129 (6). https://doi.org/10.1029/2024JF007667 e2024JF007667.

Samimi, S., Marshall, S.J., MacFerrin, M., 2020. Meltwater penetration through temperate ice layers in the percolation zone at DYE-2, Greenland ice sheet. Geophys. Res. Lett. 47. https://doi.org/10.1029/2020GL089211 e2020GL089211.

Schyberg, H., Yang, X., Køltzow, M.A.O., Amstrup, B., Bakketun, Å., Bazile, E., Bojarova, J., Box, J.E., Dahlgren, P., Hagelin, S., Homleid, M., Horániy, A., Hoyer, J., Johansson, Å., Killie, M.A., Körnich, H., Le Moigne, P., Lindskog, M., Manninen, T., Nielsen Englyst, P., Nielsen, K.P., Olsson, E., Palmason, B., Peralta Aros, C., Randriamapiania, R., Samuelsson, P., Stappers, R., Støylen, E., Thorsteinsson, S., Valkonen, T., Wang, Z.Q., 2020. Arctic regional reanalysis on single levels from 1991 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). https://doi.org/10.24381/cds.713858f6 [dataset].

Shang, X.-Y., Zheng, L., Liang, Q., Li, T., Cheng, X., 2025. Buried depth of retained meltwater in polar ice sheets detected by multi-source microwave remote sensing observations. Adv. Clim. Chang. Res. 16 (4), 813–829. https://doi.org/10.1016/j.accre.2025.05.006.

Singh, S., Durand, M., Kim, E., Barros, A.P., 2024. Bayesian physical-statistical retrieval of snow water equivalent and snow depth from X- and Ku-band synthetic aperture radar - demonstration using airborne SnowSAR in SnowEx'17. Cryosphere 18, 747–773. https://doi.org/10.5194/tc-18-747-2024.
Steffen, K., Box, J. E., and Abdalati, W., 1996. Greenland climate network: GC-Net. In: Glaciers, Ice Sheets and Volcanoes: A Tribute to Mark F. Meier, edited by: Colbeck, S., US Army Cold Regions Reattach and Engineering (CRREL), CRREL Special Report, 98–103, URL: https://apps.dtic.mil/sti/citations/ADA321342.

Steger, C.R., Reijmer, C.H., van den Broeke, M.R., 2017. The modelled liquid water balance of the Greenland ice sheet. Cryosphere 11, 2507–2526. https://doi.org/10.5194/tc-11-2507-2017.

Tedstone, A.J., Machguth, H., 2022. Increasing surface runoff from Greenland's firn areas. Nat. Clim. Chang. 12, 672–676. https://doi.org/10.1038/s41558-022-01371-z.

Tedstone, A., Machguth, H., Clerx, N., et al., 2025. Concurrent superimposed ice formation and meltwater runoff on Greenland's ice slabs. Nat. Commun. 16, 4494. https://doi.org/10.1038/s41467-025-59237-9.

Torinesi, O., Fily, M., Genthon, C., 2003. Interannual variability and trend of the Antarctic summer melting period from 20 years of spaceborne microwave data. J. Clim. 16, 1047–1060. https://doi.org/10.1175/1520-0442(2003)016<1047:VATOTS>2.0.CO;2.

Van As, D., van den Broeke, M., Reijmer, C., van de Wal, R., 2005. The summer surface energy balance of the high Antarctic plateau. Boundary-Layer Meteorol. 115, 289–317. https://doi.org/10.1007/s10546-004-4631-1.

van den Broeke, M.R., Enderlin, E.M., Howat, I.M., Kuipers Munneke, P., Noël, B.P.Y., van de Berg, W.J., van Meijgaard, E., Wouters, B., 2016. On the recent contribution of the Greenland ice sheet to sea level change. Cryosphere 10, 1933–1946. https://doi.org/10.5194/tc-10-1933-2016.

van den Broeke, M., Kuipers Munneke, P., Noël, B., Reijmer, C.H., Smeets, P., van de Berg, W.J., van Wessem, M., 2023. Contrasting current and future surface melt rates on the ice sheets of Greenland and Antarctica: lessons from in situ observations and climate models. PLOS Clim. 2, e0000203. https://doi.org/10.1371/journal.pcm.0000203.

[dataset] Vandecrux, B., 2025a. Depth of subsurface water on the Greenland ice sheet from machine learning and multifrequency passive microwave measurements (2010–2023). GEUS Dataverse. https://doi.org/10.22008/FK2/69CRRT.

Vandecrux, B., 2025b. Modeled brightness temperatures at 1.4–18.7 GHz at 19 Greenland ice sheet sites using SMRT and the GEUS snow model. GEUS Dataverse. https://doi.org/10.22008/FK2/KKIOZZ [dataset].

[software] Vandecrux, B., Picard, G., 2025. Estimating Subsurface Liquid Water Depth on the Greenland Ice Sheet Using Passive Microwave Observations, Snow Modeling, and Machine Learning (v1.0). Zenodo. https://doi.org/10.5281/zenodo.15791029.

Vandecrux, B., Fausto, R.S., Langen, P.L., van As, D., MacFerrin, M., Colgan, W.T., Ingeman-Nielsen, T., Steffen, K., Jensen, N.S., Møller, M.T., Box, J.E., 2018. Drivers of firn density on the Greenland ice sheet revealed by weather station observations and modeling. Journal of Geophysical Research: Earth Surface 123, 2563–2576. https://doi.org/10.1029/2017JF004597.

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., Box, J.E., 2019. Firn data compilation reveals widespread decrease of firn air content in western Greenland. The Cryosphere 13, 845–859. https://doi.org/10.5194/tc-13-845-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., Box, J. E., 2020a. Firn cold content evolution at nine sites on the Greenland ice sheet between 1998 and 2017. J. Glaciol. 66, 591–602. https://doi.org/10.1017/jog.2020.30.

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., Box, J.E., 2020b.

The firn meltwater retention model intercomparison project (RetMIP): evaluation of nine firn models at four weather station sites on the Greenland ice sheet. Cryosphere 14, 3785–3810. https://doi.org/10.5194/tc-14-3785-2020.

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., Steffen, K., 2023. The historical Greenland climate network (GC-net) curated and augmented level-1 dataset. Earth Syst. Sci. Data 15, 5467–5489. https://doi.org/10.5194/essd-15-5467-2023.

Vandecrux, B., Fausto, R.S., Box, J.E., Covi, F., Hock, R., Rennermalm, Å.K., Heilig, A., Abermann, J., van As, D., Bjerre, E., Fettweis, X., Smeets, P.C.J.P., Kuipers Munneke, P., van den Broeke, M.R., Brils, M., Langen, P.L., Mottram, R., Ahlstrøm, A.P., 2024a. Recent warming trends of the Greenland ice sheet documented by historical firn and ice temperature observations and machine learning. Cryosphere 18, 609–631. https://doi.org/10.5194/tc-18-609-2024.

[dataset] Vandecrux, B., Amory, C., Ahlstrøm, A.P., Akers, P.D., Albert, M., Alley, R.B., Alves de Castro, M., Arnaud, L., Baker, I., Bales, R., Benson, C., Box, J.E., Brucker, L., Buizert, C., Chandler, D.M., Charalampidis, C., Cherblanc, C., Clerx, N., Colgan, W., Covi, F., Dattler, M., Denis, G., Derksen, C., Dibb, J.E., Ding, M., Dixon, D., Eisen, O., Fahrner, D., Fausto, R., Favier, V., Fernandoy, F., Freitag, J., Gerland, S., Harper, J., Hawley, R.L., Heuer, J., Hock, R., Hou, S., How, P., Hubbard, B., Humphrey, N., Iizuka, Y., Isaksson, E., Kameda, T., Karlsson, N.B., Vladimirova, D., Kjær, H.A., Kreutz, K., Kuipers Munneke, P., Lazzara, M.A., Lemeur, E., Lenaerts, J.T.M., Lewis, G., Lindau, F.G.L., Lindsey-Clark, J., MacFerrin, M., Machguth, H., Magand, O., Mankoff, K.D., Marquetto, L., Martinerie, P., Mayewski, P.A., McConnell, J.R., Medley, B., Miège, C., Miles, K.E., Miller, H., Miller, O., Montgomery, L., Morris, E., Mosley-Thompson, E., Mulvaney, R., Niwano, M., Oerter, H., Osterberg, E., Otosaka, I., Picard, G., Polashenski, C., Reijmer, C., Rennermalm, A., Rutishauser, A., Scanlan, K.M., Simoes, J.C., Simonsen, S.B., Smeets, P.C., Smith, A., Solgaard, A., Spencer, M., Steen-Larsen, H.C., Stevens, C.M., Sugiyama, S., Svensson, J., Tedesco, M., Thomas, E.R., Thompson-Munson, M., Tsutaki, S., van As, D., Van den Broeke, M.R., van Tiggelen, M., Wang, Y., Wilhelms, F., Winstrup, M., Xiao, J., Xiao, C., 2024b. The SUMup collaborative database: surface mass balance, subsurface temperature and density measurements from the Greenland and Antarctic ice sheets (2024 release). Arctic Data Center. https://doi.org/10.18739/A2M61BR5M.

Yang, J., Jiang, L., Luojus, K., Pan, J., Lemmetyinen, J., Takala, M., Wu, S., 2020. Snow depth estimation and historical data reconstruction over China based on a random forest machine learning approach. Cryosphere 14, 1763–1778. https://doi.org/10.5194/tc-14-1763-2020.

Zeiger, P., Picard, G., 2024a. SMOS rSIR-enhanced brightness temperature maps on Greenland and Antarctica (v1.0). EaSy Data. https://doi.org/10.57932/t72f9515-3699-4fae-92fc-a350075d042f [dataset].

Zeiger, P., Picard, G., 2024b. SMOS rSIR-enhanced melt maps on Greenland and Antarctica (v1.0). EaSy Data. https://doi.org/10.57932/1970fb7c-cdb4-4ddf-9891-1-e6836a46f25 [dataset].

Zeiger, P., Picard, G., Richaume, P., Mialon, A., Rodriguez-Fernandez, N., 2024. Resolution enhancement of SMOS brightness temperatures: application to melt detection on the Antarctic and Greenland ice sheets. Remote Sens. Environ. 315, 114469. https://doi.org/10.1016/j.rse.2024.114469.

Zheng, L., Cheng, X., Shang, X., Chen, Z., Liang, Q., Wang, K., 2022. Greenland ice sheet daily surface melt flux observed from space. Geophys. Res. Lett. 49. https://doi.org/10.1029/2021GL096690 e2021GL096690.

Zwally, J.H., 1977. Microwave emissivity and accumulation rate of polar firn. J. Glaciol. 18, 195–215. https://doi.org/10.3189/S0022143000021304.

Zwally, J.H., Fiegles, S., 1994. Extent and duration of Antarctic surface melting. J. Glaciol. 40, 463–475. https://doi.org/10.3189/S0022143000012338.