Runoff from Greenland's firn area – why do MODIS, RCMs and a firn model disagree?

Due to increasing air temperatures, surface melt and meltwater runoff expand to ever higher elevations on the Greenland ice sheet and reach far into its firn area. Here, we evaluate how two regional climate models (RCMs) simulate the expansion of the ice sheet runoff area: MAR, and RACMO with its offline firn model IMAU-FDM. For the purpose of this comparison we first improve an existing algorithm to detect daily visible runoff limits from MODIS satellite imagery. We then apply the improved algorithm to most of the Greenland ice sheet and compare MODIS to RCM runoff limits for the years 2000 to 2021. We find that RACMO/IMAU-FDM runoff limits are on average somewhat lower than MODIS and show little fluctuation from year to year. MAR runoff limits are higher than MODIS, but their inter-annual fluctuations are more similar to MODIS. Both models apply a bucket scheme to route meltwater vertically. Focusing on the K-transect, we demonstrate that differences in modelled firn temperatures and in the implementation of the bucket scheme govern RCM simulated runoff limits. The formulation of the runoff condition is of large influence: in RACMO/IMAU-FDM meltwater is only considered runoff when it reaches the bottom of the simulated firn pack; in MAR runoff can also occur from within the firn pack, which contributes to its high runoff limits. We show that total runoff along the K-transect, simulated by the two RCMs, diverges by up to 29 % in extraordinary melt years. This difference is mostly caused by the diverging simulated runoff limits, which emphasizes the importance of improving the simulations of Greenland's melting firn area.

1 Introduction

Polar regional climate models (RCMs) are widely used to assess past, present and future surface mass balance of the Greenland and Antarctic ice sheets (Box et al., 2004; Fettweis et al., 2008; Noël et al., 2016; IMBIE Team, 2018, 2020). The accuracy of RCM output relies, among other factors, on data available for model calibration and evaluation. Essential for RCM evaluation are meteorological observations (e.g. Steffen and Box, 2001; Fausto et al., 2021), surface mass balance measurements (e.g. Benson, 1962; Greuell et al., 2001; van de Berg et al., 2006; Machguth et al., 2016b; Karlsson et al., 2016; Fausto et al., 2021) and remote sensing products (e.g. Fettweis et al., 2006; Mohajerani et al., 2019; Slater et al., 2021). RCMs have been extensively evaluated for Greenland's ablation area (e.g. Gallee and Duynkerke, 1997; Lefebre et al., 2005; Noël et al., 2016) and its higher accumulation area (e.g. Rae et al., 2012; Noël et al., 2016) and have been found to perform well when compared to meteorological observations, surface mass balance measured at stake locations as well as in ice cores and gravimetric ice sheet mass balance (e.g. Fettweis et al., 2017, 2020).

Comprehensive model evaluation requires also testing the RCMs in the transition zone in-between the ablation and the higher accumulation area. In this area a delicate balance exists between accumulation and ablation processes. In summer, when melt, runoff and accumulation can occur simultaneously, working conditions are challenging (e.g. Holmes, 1955; Clerx et al., 2022). Consequently, the availability of field data is limited and few studies (e.g. Covi et al., 2022; Zhang et al., 2023; Vandecrux et al., 2024) have evaluated RCMs in this transition zone.

Within the elevation range of the transition zone lie the equilibrium line altitude (ELA) and the runoff limit. The former is the elevation where the climatic mass balance equals zero (Cogley et al., 2011), the latter is defined as the uppermost elevation from where meltwater can reach the ocean and contribute to mass loss (Cogley et al., 2011; Tedstone and Machguth, 2022; Clerx et al., 2022). The elevation of the equilibrium line and the runoff limit varies from year to year in response to weather conditions. Thereby, the runoff limit lies within the accumulation area (Shumskii, 1955, 1964) and is thus located above the ELA.

Tedstone and Machguth (2022) compared seasonal maxima of visible runoff limits mapped from Landsat satellite imagery to runoff extent simulated by the two RCMs RACMO 2.3p2 (Noël et al., 2018) and MAR v3.11 (Fettweis et al., 2017) forced by ERA-40/ERA-I/ERA5. The comparison revealed substantial differences between RCMs and remotely sensed visible runoff limits, but also between the two RCMs involved. While remotely sensed visible runoff limits are subject to uncertainties, it remains unclear what causes the remarkable differences between the RCMs. If RCMs differ in simulating the runoff extent of the Greenland ice sheet, this results in inaccuracies in future scenarios of mass loss and sea-level contribution. Indeed, Glaude et al. (2024) found large differences in RCM simulated runoff area for the year 2100 under a high-end warming scenario (SSP5-8.5). Glaude et al. (2024) point out that the three RCMs studied, among them RACMO and MAR, differ by a factor of two in their predicted surface mass balance for the year 2100.

Here we aim at explaining why simulated runoff limits differ between models. For this purpose we compare remotely sensed visible runoff limits and simulated runoff limits by MAR, RACMO and the firn model IMAU-FDM. We use daily Moderate Resolution Imaging Spectroradiometer (MODIS) visible runoff limits for the years 2000 to 2021, derived by an improved version of the algorithm used in Machguth et al. (2023). We use MODIS visible runoff limits instead of the aforementioned Landsat visible runoff limits because MODIS offers higher temporal resolution. We analyze the differences between observed and modelled runoff limits in the context of modelled parameters that potentially influence simulated runoff. Among the selected parameters are surface albedo, firn density and temperature, as well as refreezing. We identify which of the parameterizations in the models likely cause the deviations. Finally, we quantify their impact on simulated mass balance along a transect in south-west Greenland.

2 Data

2.1 Data for MODIS visible runoff limit detection

The detection of MODIS visible runoff limits Υ_obs is based on an optimized version of the algorithm by Machguth et al. (2023). The improved algorithm (Sect. 3.1) relies on the following input: (i) daily MOD10A1 data (MODIS/Terra Daily Snow Cover at 500 m resolution, version 6.0; Hall and Riggs, 2016); (ii) daily MOD09GA data (MODIS/Terra Surface Reflectance Daily at 500 m, version 6.0; Vermote and Wolfe, 2015); (iii) the Arctic DEM (100 m resolution mosaic, v.3.0; Porter et al., 2018, here downsampled to the 500 m MODIS grid); (iv) outlines of the Greenland ice sheet according to Rastner et al. (2012) and (v) Greenland-wide arrays of surface ice flow velocity in x and y direction (Joughin et al., 2016, 2017).

2.2 Model data

To quantify modelled runoff limits Υ_rcm we use (i) simulated runoff from the polar regional climate model MAR (version 3.14, 10 km resolution, forced by ERA5), (ii) the polar version of the Regional Atmospheric Climate Model RACMO (Noël et al., 2019, version 2.3p2 on the grid FGRN055, forced by ERA5) at a resolution of 5.5 km as well as (iii) the offline firn model IMAU-FDM v1.2G (Ligtenberg et al., 2011, 2018; Brils et al., 2022). Descriptions of these three models, with special focus on their firn simulation, are provided in Sect. 3.2. RACMO data are frequently used in a version that is further downscaled to 1 km resolution and bias corrected (Noël et al., 2016, 2019). The downscaled data have a temporal resolution of 1 d, which is insufficient to force IMAU-FDM v1.2G, so for firn applications these data cannot be used.

We use a set of RCM parameters (Table 1) to explore the reasons behind potential differences in MODIS and RCM runoff limits. Various parameters are not written to output by RACMO2.3p2 and are thus unavailable. Instead, we obtained them from the offline firn model IMAU-FDM v1.2G henceforth IMAU-FDM. The model is forced in offline mode by RACMO2.3p2 and is run on an identical spatial grid. In the following our usage of the term RCM also refers to the offline firn model IMAU-FDM. As explained in Sect. 3.2.2, the latter is very similar to RACMO's firn module. Furthermore, we refer to “MAR” for MARv3.14 and to “RACMO” for RACMO2.3p2.

Table 1 List of RCM simulated parameters used to calculate and investigate runoff limits. “RACMO” stands for RACMO2.3p2 at native 5.5 km resolution, “IMAU-FDM” stands for IMAU-FDM v1.2G and “MAR” stands for MARv3.14.

MAR output is obtained at daily temporal resolution. Output from RACMO and IMAU-FDM are at 10 d intervals. Where needed, MAR data are averaged or summed to the lower temporal resolution.

3 Methods

3.1 Detecting MODIS Υ_obs along flowlines

The algorithms by Greuell and Knap (2000) and Machguth et al. (2023) detect Υ_obs on AVHRR (1.1 km spatial resolution; Greuell and Knap, 2000) or MODIS (500 m resolution; Machguth et al., 2023) satellite imagery. Given the low spatial resolution as compared to e.g. Landsat, Υ_obs is identified indirectly, that is where spatial variability of surface albedo α transitions from low to high. Low spatial variability of α indicates a monotonous snow covered surface. Variability of α is high where dark meltwater streams, lakes and slush fields intersect the bright snow cover. Despite this indirect approach, MODIS Υ_obs highly agree with visible runoff limits detected on finer resolution (30 m) Landsat imagery (Machguth et al., 2023).

Machguth et al. (2023) scanned rectangular polygons of width p_w and length p_l≫p_w for the location where the standard deviation of surface albedo σ_α falls below a certain threshold. If a set of additional conditions and tests are fulfilled (see Machguth et al., 2023), the location is considered to represent Υ_obs. The long axes of the polygons needed to be oriented along the strongest gradient in α, which is in the direction of the surface slope. Polygons in Machguth et al. (2023) were strictly oriented west-east. Consequently, the application of the method was restricted to areas of the western flank of the ice sheet.

Here we apply the method by Machguth et al. (2023) with two major modifications that allow application to all of the Greenland Ice Sheet: (1) We create so called flowline-polygons of p_w=20 km, henceforth simply called flowlines, and (2) implement an improved calculation of σ_α. The former allows detection of Υ_obs in complex topography sloping in any direction, the latter improves detection of Υ_obs by calculating and subtracting the influence of temporally persistent albedo features. These modifications, as well as smaller optimizations, are detailed in Appendix A and Fig. A1. For further details on the algorithm we refer to Machguth et al. (2023).

3.2 RCM simulations of the firn cover and runoff

3.2.1 MAR

We use daily outputs at 10 km resolution from version 3.14 of MAR, forced every 6 h by the ERA5 reanalysis. The data are composed of two transient simulations: the first one starts in September 1974 but only the period 1980–1999 is used. The second one begins in September 1994 and the period 2000–2023 is used. Together, the two simulations cover the years 1980 to 2023. In the set-up used here, MAR resolves the uppermost 21 m of snow and firn using a time-varying number of layers up to a maximum of 21 layers. For densities lower than 450 kg m⁻³, the CROCUS snow model albedo (Brun et al., 1992) is used with a minimum value of 0.7. Where surface density exceeds 450 kg m⁻³, the minimum value of albedo declines between the minimum snow albedo (0.7) and clean ice albedo (0.55) as a linear function of increasing density. On bare ice (surface density higher than 900 kg m⁻³), CROCUS snow model albedo is not used and the albedo varies exponentially between 0.55 (clean ice) and 0.5 (wet ice) as a function of the accumulated surface water height and the slope. The dependency on water depth and slope follows Lefebre et al. (2003) but the albedo for large water depths is set to 0.5 instead of the original 0.15. The parametrization is of limited impact but is maintained to address the effect of supraglacial lakes in future model versions.

The main changes of MARv3.14 with respect to MARv3.12 (Vandecrux et al., 2024) are as follows: Some bugs in the clouds scheme have been corrected and a continuous snowfall-rainfall limit has been introduced for near-surface temperature between −1 °C (100 % of precipitation falls as snow) and +1 °C (100 % rain). MARv3.14 now uses the radiative scheme from ERA5 (Hogan and Bozzo, 2018; Grailet et al., 2025) instead of the one from ERA40 (Morcrette, 2002) in former MAR versions.

MAR parameterizes meltwater percolation through an instantaneous bucket scheme. Slush is not allowed in MARv3.14 simulations and the maximum liquid water saturation in snow and firn (i.e. irreducible water saturation, expressed in % of the pore volume) is 7 % at the surface and linearly reduces to 2 % at 1 m depth. Below that depth, irreducible water saturation is set to 2 % (earlier MAR versions used higher irreducible water saturations; see e.g. Alexander et al., 2019). Meltwater that percolates into a snow or firn layer can refreeze if the layer temperature is below 0 °C or it can be retained as irreducible water if the layer is temperate. If neither of the two processes are possible, that is if the layer has become temperate and irreducible water saturation is at its maxima, the remaining meltwater will either percolate to the next layer below or run off immediately. The following conditions decide between percolation and immediate runoff. If the density of a layer is <830 kg m⁻³, percolation to the next deeper layer takes place. For layers of density ≥830 kg m⁻³, a density runoff threshold determines how much of any meltwater gets removed immediately as runoff: 0 % for 830 kg m⁻³ to 100 % for densities above 900 kg m⁻³. The remainder percolates to the next layer below. Where ice lenses are simulated by MAR, $\mathrm{2}/\mathrm{3}$ of the percolating meltwater progress to underlying layers and the remaining $\mathrm{1}/\mathrm{3}$ are considered run off. Thereby an ice lens is defined as a layer with a density of >900 kg m⁻³ that lies on top of a layer where density is ≤900 kg m⁻³. Furthermore, any meltwater that reaches the bottom of the MAR firn column is also considered runoff.

For further details on MAR we refer to Fettweis et al. (2013, 2017, 2020). Previous MAR versions have been successfully validated over the Greenland ice sheet by comparison with surface mass balance measurements (Fettweis et al., 2020), satellite derived melt extent (Fettweis et al., 2011) and in situ atmospheric measurements (Delhasse et al., 2020).

3.2.2 IMAU-FDM and RACMO

We use data from the RCM RACMO and the offline firn model IMAU-FDM, which is very similar to RACMO's firn module. While it would be preferable to consistently use RACMO data for comparison to MAR, we here use IMAU-FDM firn simulations because RACMO outputs only depth integrated firn data. In the following we first explain IMAU-FDM, then explain differences to RACMO's firn module, and finally provide information on RACMO and its forcing of IMAU-FDM.

IMAU-FDM v1.2G (Brils et al., 2022) is a semi-empirical firn densification model that simulates the time evolution of firn density, temperature, liquid water saturation and changes in surface elevation owing to variability of firn depth. Vertical water transport in IMAU-FDM is instantaneous and calculated via the bucket method. When liquid water is added to the firn column by melt or rain, it is transported vertically downwards. Starting at the uppermost model layer, the scheme checks if there is cold content and pore space available for refreezing. If so, refreezing takes place, raising the layer's temperature and density, until either (i) all water has been refrozen, (ii) the layer has turned into ice (i.e. has a density of 917 kg m⁻³), or (iii) reaches 0 °C. Irreducible water will be retained in liquid form within the pores of a temperate firn layer. The maximum amount that can be retained depends on the layer's porosity, following Coléou and Lesaffre (1998) (irreducible water saturation is ∼5.8 % at a snow density of 300 kg m⁻³ and ∼15 % at 800 kg m⁻³. Any water that cannot refreeze or be retained as irreducible water will percolate to the next layer below. These steps are then repeated in the next firn layer, and so on until no more liquid water is present aside the irreducible water saturation within temperate layers. This is all done within a single time step, which means that vertical percolation is instantaneous. The bucket method also implies that liquid water percolates through any ice layer, because they contain no pore space to accommodate refreezing. Water is not allowed to pond or run off on top of ice layers.

When the water reaches the interface between firn and glacial ice, it is assumed to run off instantaneously. The depth of the horizontal modelling domain of IMAU-FDM varies in space and time and is defined by the condition that the deepest 200 grid cells must all exceed a density of 910 kg m⁻³. Consequently the thickness of the firn layer, that is from the surface to the depth below which all grid cells exceed a density of 830 kg m⁻³, varies and reaches maxima of 100 m in high-accumulation regions of the south-east of the ice sheet. A more typical maximum firn thickness is ∼70 m.

RACMO's firn module also simulates the firn column from the surface down to glacial ice and uses similar physical parametrisations as IMAU-FDM, albeit at a lower vertical resolution (max. 100 but typically 40 layers in RACMO; up to 3000 layers in IMAU-FDM) and less comprehensive initialisation to save computing costs.

IMAU-FDM is forced at the upper boundary by 3-hourly RACMO surface temperature and mass fluxes, interpolated to 15 min. In RACMO, the snow albedo scheme is based on prognostic snow grain size, cloud optical thickness, solar zenith angle and impurity concentration in snow (Kuipers Munneke et al., 2011). Impurity concentration is assumed constant in time and space. Bare ice albedo is prescribed from the 500 m MODIS 16 d albedo version 5 product (MCD43A3v5) as the lowest 5 % surface albedo records for the period 2000–2015. Thresholds are applied to these values: minimum ice albedo is set to 0.3 for dark ice in the low-lying ablation zone, and a maximum value of 0.55 is used for bright ice under perennial snow cover in the accumulation zone, i.e. only used when all firn melts away which does not happen in this run. RACMO snow albedo typically ranges between ∼0.7 for highly metamorphosed, coarse grained snow under clear-sky conditions and ∼0.95 for fine grained snow under cloudy conditions. RACMO2.3p2 surface energy balance, surface mass balance and melt output over the GrIS have been extensively evaluated, notably along the K-transect, and were found to be generally robust (Noël et al., 2019).

3.3 Calculating Υ_rcm from RCM output

We distinguish between daily runoff limits Υ_rcm and annual maximum runoff limits maxΥ_rcm, which mark the highest elevation where runoff occurs for each year. Both Υ_rcm and maxΥ_rcm are calculated on the same 20 km wide flowlines as used for the detection of Υ_obs. For each flowline, we consider RCM grid cells whose center falls within the flowline. Given the elevation of each grid cell and simulated runoff, we then calculate runoff against elevation.

There is no generally accepted definition of maxΥ_rcm in terms of runoff per year. Tedstone and Machguth (2022) quantified the sensitivity of maxΥ_rcm to runoff thresholds of >1, >5, >10, and >20 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$. They found that MAR and RACMO maxΥ_rcm are rather insensitive to the choice of threshold. Furthermore, they stated that the uncertainties associated with the choice of thresholds are small compared to the substantial differences in maxΥ_rcm between the two RCMs. We here adopted their chosen threshold of >10 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{a}}^{-\mathrm{1}}$ to calculate maxΥ_rcm. To estimate daily Υ_rcm we use a threshold of >1 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{d}}^{-\mathrm{1}}$.

3.4 Analyzing RCM process simulations near the runoff limit

Our goal is to understand why deviations occur (i) between Υ_obs and Υ_rcm and (ii) between the two Υ_rcm (labeled ${\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ and ${\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$). We focus this part of the analysis on the K-transect which has been studied intensively with respect to ice sheet boundary layer meteorology (van den Broeke et al., 1994), surface mass balance (Van de Wal et al., 2005, 2012), firn processes (Machguth et al., 2016a; Mikkelsen et al., 2016; Rennermalm et al., 2021) and firn hydrology (Clerx et al., 2022). Here we defined the K-transect as the line that follows the 67° N parallel, starts at the ice margin at ∼250 m a.s.l. / 50° W, and reaches to the ice divide at ∼2520 m a.s.l. / 42.7° W (Fig. 1). For both RCMs and IMAU-FDM, we extract the grid cells which are closest to the ∼320 km long transect. This results in lines of RCM grid cells which are one cell wide and 33 (MAR) or 57 cells (RACMO, IMAU-FDM) in length.

Figure 1 Median, highest and lowest of all annual MODIS maxΥ_obs for the time period 2000 to 2021. Retrievals at the east coast between 60 to 68.4° N, where the hydrological regime is dominated by firn aquifers, have been masked. Flowlines highlighted in orange indicate the locations for which detailed results are shown in Figs. 2, 4, C1 and C2. The location of the K-transect (Figs. 5, 6 and C3) is indicated as well.

Along the K-transect we analyse the RCM simulated parameters listed in Table 1. We graphically display temporal and spatial changes and visually search for parameters that show peculiar or unexpected values in the broader elevation range around the runoff limit. If found, we investigate the underlying RCM parameterizations in order to understand their potential influence on Υ_rcm.

4 Results

4.1 MODIS Υ_obs detections

Figure 1 summarizes the MODIS-derived Υ_obs for all of the Greenland Ice Sheet. The approach creates few and meaningless Υ_obs in areas dominated by meltwater discharge through aquifers. This is to be expected as surface meltwater features are largely absent in such areas. Consequently Fig. 1 does not show retrievals from 60 to 68.4° N along Greenland's east coast. However, we show detected Υ_obs located in smaller aquifer regions elsewhere on the ice sheet. Excluding retrievals from 60 to 68.4° N along the east cost, 63 400 Υ_obs in 417 flowlines remain, which corresponds on average to ∼7 retrievals per flowline and year. The actual number of annual retrievals varies geographically and is highest in the southwest, exceeding on average 18 retrievals per flowline and melt season.

Compared to Machguth et al. (2023) and their study area, we find that the updated algorithm yields ∼80 % more Υ_obs detections. This difference is mainly due to the new algorithm being able to place more flowlines that are optimized for complex topographies. The average number of Υ_obs detections per flowline is 5.5 % higher than per stripe, which were the strictly east-west oriented bands in Machguth et al. (2023). Outside of the area investigated by Machguth et al. (2023), the new approach provides numerous detections of Υ_obs in the north-west of Greenland, from near Pituffik Space Base to Humboldt and Petermann glaciers, as well as in the region of the north-east Greenland ice stream. Few detections occur along the central part of the east coast where the terrain is complex and steep, with numerous outlet glaciers. The approach appears not well suited to such terrain because most outlet glaciers are narrow, compared to the 20 km width of the flowline polygons. Consequently, along the outlet glaciers few glacier pixels are available for retrieval of the Υ_obs. Apart from Petermann Glacier, there are few detections beyond 80° N, the reasons for which are unclear. Tedstone and Machguth (2022), who used Landsat to detect surface hydrology, also noted few detections in the region.

Figure A2 compares Υ_obs to the Landsat-derived visible runoff limits from Tedstone and Machguth (2022). The comparison yields a good agreement between the two data sets and is discussed in Appendix A4.

Figures 2, C1 and C2 exemplify the temporal detail of the Υ_obs data. The figures demonstrate frequent behavior where Υ_obs rises relatively early in the melt season and reaches a plateau before melting ends (see also Machguth et al., 2023). By design of the detection and filtering algorithms, there is typically no decrease in Υ_obs towards the end of the melt season: Most decreasing Υ_obs are filtered out because optical remote sensing is poorly suited to detect continued hydrological activity under freshly fallen autumn snow (Machguth et al., 2023).

Figure 2 Evolution of MODIS visible runoff limits Υ_obs and RCM simulate runoff limit Υ_rcm over two selected melt seasons and flowlines. Solid lines show RCM runoff limits at daily resolution for MAR and in 10 d steps for IMAU-FDM. Subplot (a) shows data for the transect NE for the year 2010, (b) shows the transect CW for the year 2009. See Fig. 1 for the location of the two transects shown.

4.2 Comparing Υ_obs and Υ_rcm

4.2.1 Comparing annual maxima

Figure 3 shows how maxΥ_obs and maxΥ_rcm vary along Greenland's western flank. The RCMs and MODIS show a general decrease of the runoff limit towards higher latitudes (Fig. 3b). Certain deviations from this trend are common to all data: maxΥ_obs and maxΥ_rcm are depressed south of ∼63° N and elevated in-between ∼71 and ∼72.5° N. Where firn aquifers are present, maxΥ_obs are biased low and standard deviation is increased. Otherwise, the differences between maxΥ_obs and maxΥ_rcm depend strongly on the RCM. IMAU-FDM simulated runoff limits are on average at 1545 m a.s.l., lower than maxΥ_obs which are on average at 1613 m a.s.l. IMAU-FDM have a low standard deviation of 31 m compared to MODIS (99 m). MAR maxΥ_rcm are at 1816 ± 94 m a.s.l., substantially higher than MODIS but with similar standard deviation. Figure 4 illustrates for two selected regions how maxΥ_obs and maxΥ_rcm fluctuate over time. IMAU-FDM simulated runoff limits vary little between the years. The intense melt seasons of 2012 and 2019 leave virtually no trace in its runoff limits. MAR maxΥ_rcm vary with the intensity of the melt season. Temporal variability of $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ exceeds MODIS in the south (Fig. 4b), but is rather similar further north (Fig. 4a).

Figure 3 The western slope of the Greenland ice sheet and mean MODIS, MAR and IMAU-FDM runoff limits, averaged over the time period 2000 to 2021. (a) Map of Greenland's west coast showing the flowlines along which the runoff limits have been calculated. (b) Mean and standard deviation of MODIS maxΥ_obs and maxΥ_rcm of MAR and IMAU-FDM for all flowlines that fall into the area shown. Gray shading indicates latitudes where firn aquifers occur (Miège et al., 2016).

Figure 4 Annual mean of MODIS maxΥ_obs and maxΥ_rcm of MAR and IMAU-FDM. (a) Averaged over the six flowlines of region NW and (b) averaged over the six flowlines at around the K-transect (region SW, see Fig. 1). Shading illustrates annual variability (±1σ) of maxΥ_obs or maxΥ_rcm within the two groups of six neighboring flowlines.

4.2.2 Comparing seasonal evolution of Υ_obs and Υ_rcm

Comparing the seasonal evolution of Υ_rcm and Υ_obs shows that MODIS and RCM runoff limits often reach their seasonal maxima at similar points in time (Figs. 2, C1 and C2). The dates of the first appearance of the runoff limit are often similar between RCMs and MODIS. However, Υ_rcm fluctuate strongly, often dropping and increasing, within a few days, over hundreds of meters in elevation (e.g. Fig. 2). The effect is more pronounced for MAR which is due to the higher temporal resolution of the MAR data. MODIS Υ_obs indicate a more continuous process where the visible runoff limit remains at high elevations, also during cold spells.

Agreement of ${\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ to the seasonal evolution of Υ_obs is generally good (Figs. 2, C1 and C2). However, ${\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ always tends to reach its maxima at very similar elevations, regardless of the intensity of the melt season. This is the same behavior shown for $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ in Fig. 4. MAR Υ_rcm typically overshoot Υ_obs (Figs. C1 and C2).

4.3 RCM process simulations at the runoff limit

Potential causes for the large differences between Υ_rcm are (i) differences in the amount of simulated melt or snowfall in MAR or RACMO, or (ii) differences in the firn parameterizations that impact simulated runoff. In Appendix B we demonstrate that differences in melt or accumulation at the maxΥ_rcm are small and cannot explain the differences between $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ and $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$. Here we therefore investigate whether reasons for the differences in maxΥ_rcm can be found in the models' firn parameterizations. For the sake of clarity, we focus the analysis on the K-Transect, whose representativeness for the entire ice sheet will be assessed in the Discussion. Furthermore, we focus on the two contrasting melt seasons of 2012 and 2017. The former was dominated by early, persistent and intense melting, the latter by intermittent and moderate melt. They represent the end members of the last 25 mass balance years that were dominated by mass loss (see Fig. B1).

Figures 5 and 6 visualize and compare RCM simulated parameters for the 2012 and 2017 melt seasons. Figure 5 shows average or summed values over the time period 1 May to 31 October and Fig. 6 illustrates the spatio-temporal evolution of parameters over the same time frame. In 2012, IMAU-FDM shows discontinuities at the location of $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$: Mean albedo increases by ∼0.05 (Fig. 5c) while melt drops by ∼400 mm w.e. or 31 % (Fig. 5e). The contrast in albedo is even higher (an increase from 0.65 to 0.78) when averaging only from mid-July to mid-August 2012. At $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$, runoff drops from slightly higher than 1000 mm w.e. to zero (Fig. 5e). Across $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$, the percentage of melt running off drops from ∼80 % to zero (Fig. 5g). This sudden shut-down of runoff is compensated by an abrupt increase in refreezing (Fig. 5i). In 2012 these transitions take place over the distance of a single grid cell (5.5 km), whereas in 2017, IMAU-FDM shows gradual transitions without discontinuities. In 2012, MAR shows no discontinuities in albedo and melt across $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ (Fig. 5c and e) but it exhibits step-wise changes in runoff and refreezing (Fig. 5g and i). These discontinuities are somewhat less pronounced than for IMAU-FDM. In 2017, simulated refreezing of MAR and IMAU-FDM are rather similar along the transect (Fig. 5k), regardless of $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ being located at higher elevation.

Figure 5 Comparison of RCM simulated parameters along the K-transect. The left column of subplots refers to the 2012 melt season; 2017 is to the right. The parameters shown in each row of subplots are explained in the plot titles to the left. Summed values in subplots (e) to (k) are summed over the time frames indicated at the top; runoff and refreezing are furthermore depth integrated over the first 20 m of the firn pack.

Figure 6 Comparison of IMAU-FDM (subplots c to h) and MAR (i to o) simulated parameters along the K-transect. Subplots to the left refer to the 2012 melt season; 2017 is to the right. Refreezing and liquid water content (subplots c to f and i to m) are depth integrated over the top 20 m of the firn column. Blue and pink dots denote RACMO and MAR simulated seasonal evolution of the runoff limit, respectively. Orange circles show MODIS-mapped seasonal evolution of the visible runoff limit. All heat maps are given at 10 d temporal resolution. The parameters shown in each row of subplots are explained in the plot titles to the left.

In 2012, ${\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ remained stable over an extended time period (e.g. Fig. 6c). The sharp increase in total refreezing, observed in Fig. 5i, is the result of intense refreezing that took place during the prolonged time period when the IMAU-FDM runoff limit was at its maximum (Fig. 6e). The refreezing raised 10 m firn temperatures to 0 °C (Fig. 6g), which is unique for the decade 2010 to 2020 (Fig. C3). In 2012, MAR refreezing was also focused to directly above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ (Fig. 6l), but not as clearly as IMAU-FDM. The peak in MAR summed refreezing is thus less pronounced (Fig. 5i). We notice that MAR refreezing fluctuates somewhat randomly along the transect. These fluctuations can be observed in both years and occur mainly in-between the maxΥ_rcm of the two RCMs (Fig. 5i and k). The fluctuations can also be seen in Fig. 6l and m.

In MAR, there is less influence of refreezing on 10 m firn temperatures (Fig. 6n and o) and firn temperatures below the 2012 $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ were already very close to 0 °C before the melting started. The 2012 refreezing results in moderate firn warming above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ which then persists (Fig. C3).

Figure 7 serves to assess whether maxΥ_rcm are related to simulated firn structure. In 2012, $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ coincides with the uppermost grid cell where the top 20 m of the firn consist of ice. MAR maxΥ_rcm is underlain by less dense firn and is located much higher than the uppermost grid cell of uniform ice. Furthermore, we notice that the IMAU-FDM firn profile shows an ice slab, a zone of icy firn in the top ∼5 m of the firn profile overlying material of lower density. The slab is most pronounced directly uphill of the 2012 $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$. The MAR firn profile shows a more weakly developed zone of increased near-surface density around and above the 2012 $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$.

Figure 7 Comparison of RCM simulated firn structure along the K-transect and for the year 2012. Dotted areas signify depth intervals where ρ>830 kg m⁻³ and exceeds pore close-of density. Runoff limits are also shown for the 2012 melt season.

Firn properties simulated by MAR and IMAU-FDM differ in the vicinity of the maxΥ_rcm (Fig. 7), which mandates a more detailed comparison of firn properties. Along the K-Transect, KAN_U is the optimal site for such a comparison because (i) the site is located at 1840 m a.s.l. which places it above the highest $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ and close to the average $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$, and (ii) the site features repeated measurements of firn density (Rennermalm et al., 2021) and firn temperatures (e.g. Charalampidis et al., 2016; Vandecrux et al., 2024). Figure C4 visualizes simulated MAR and IMAU-FDM firn density evolution for the top 20 m over the time period 1980–2020 at KAN_U and Fig. C5 shows simulated firn temperature profiles and a comparison to measured 10 m depth firn temperatures. IMAU-FDM firn density evolution shows annual layers getting buried and an ice slab forming in summer 2012. Afterwards, the slab gets buried under accumulating snow and firn. In contrast to this, the observed depth of the top of the ice slab (Fig. C4a) remains close to the surface. The coarser vertical resolution of the MAR outputs makes it more difficult to follow horizons as they get buried. Simulated temperatures vary strongly at the site, being close to −15 °C in IMAU-FDM and around 0 °C in MAR. The former matches measured 10 m temperatures (around −11 °C; Fig. C5c) more closely.

4.4 Υ_rcm and its relevance for RCM simulated runoff

Along the K-Transect, but also for most other regions of the ice sheet (e.g. Fig. 3), the MAR runoff zone is larger than for IMAU-FDM. The question arises to what degree this is relevant to overall runoff. On the example of the K-transect we quantify by how much total simulated runoff is influenced by $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ being at higher elevations than $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$.

For each year from 1980 to 2020 we calculate total annual RCM runoff (i) below and (ii) above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ along the K-transect (see the inset in Fig. 8). We assume unit width for the transect so runoff has the unit m³ yr⁻¹. The first value to be calculated, termed ${\int }_{\Downarrow }Q$, can be derived for both RCMs. The second value, ${\int }_{\Uparrow }Q$, can only be calculated from MAR whose maxΥ_rcm is always higher than IMAU-FDM along the K-transect.

Figure 8 Regression of 1980 to 2020 MAR simulated runoff below and above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$. Every point corresponds to one year and the two runoff values for each year are integrated along parts of the K-transect as illustrated in the inset.

Exponential regression of the two parameters ${\int }_{\Downarrow }{Q}_{\text{MAR}}$ and ${\int }_{\Uparrow }{Q}_{\text{MAR}}$ yields R²=0.83 (Fig. 8), which means the amount of MAR runoff above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ increases exponentially as a function of the MAR runoff below $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$. If MAR and IMAU total runoff below $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ limit were similar (see the following paragraph), this implies that the difference in simulated runoff between MAR and IMAU-FDM increases in high-melt seasons. The reason for the disproportional growth is that the more intense the melt season, the further apart the two maxΥ_rcm. If ${\int }_{\Uparrow }{Q}_{\text{MAR}}$ is expressed as a percentage of ${\int }_{\Downarrow }{Q}_{\text{MAR}}$, we find that for 2012 ${\int }_{\Uparrow }{Q}_{\text{MAR}}$ corresponds to 20 % of ${\int }_{\Downarrow }{Q}_{\text{MAR}}$. For the year 2017, the percentage is 3.2 % which is somewhat lower than the mean of all years (5.7 %).

For sake of clarity, the above statistics were based on Q_MAR alone. However, simulated runoff of the two RCMs below $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ are not identical. We label the area below $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ as the “common runoff area” because $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ is always situated at lower elevation than $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$. We find that over this area, MAR simulates 8370 ± 11 650 m³ yr⁻¹ (mean ± 1 std. dev.) more runoff than IMAU-FDM. Expressed in percent, MAR simulates 5.6 % ± 7.8 % more runoff over the common runoff area than IMAU-FDM. On average, the differences in RCM runoff caused by the diverging maxΥ_rcm (9430 ± 8970 m³ yr⁻¹) are similar to the differences in runoff over the common runoff area. In the extraordinary melt season of 2012, however, the influence of the differing maxΥ_rcm clearly exceeds the differences in RCM runoff over the common runoff area. In 2012, total MAR runoff along the K-transect exceeds IMAU-FDM by 72 070 m³. This corresponds to 29 % of the total 2012 IMAU-FDM runoff. Out of the total difference, 53 370 m³ or three quarters are due to MAR runoff above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$. To examine whether this finding also holds in other extreme melt seasons, we examined the runoff in 2019, which was another extraordinary melt year. Consistent with the 2012 results, in 2019, the difference in total runoff is 34 250 m³ (16 % of the 2019 IMAU-FDM runoff) out of which 26 940 m³ or almost four fifths originate from MAR runoff above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$.

5 Discussion

There are fundamental differences between runoff processes detected from remote sensing and their simulation. Optical satellite imagery primarily detects lateral runoff, visible in slush fields and meltwater streams at the surface; sub-surface runoff cannot be sensed. In contrast, current state-of-the-art firn models or RCM firn modules simulate runoff through vertical percolation alone; lateral flow is not simulated. Nevertheless, we here compared modelled and remotely sensed runoff limits on the Greenland Ice Sheet because (i) modelled runoff has the purpose of mimicking the actual, strongly lateral, process. Thus we here tested whether the mimicking approximates the effects of the actual hydrological processes. (ii) The remotely sensed visible runoff limit approximates the actual (invisible) runoff limit reasonably well at the peak of the melt season (Holmes, 1955; Clerx et al., 2022; Tedstone and Machguth, 2022).

5.1 Comparing MODIS and simulated runoff limits

We observe a relationship between maxΥ_obs and maxΥ_rcm that is in broad agreement to Ryan et al. (2019) who compared snow lines simulated by MAR, RACMO and observed from remote sensing (cf. Fig. 4 herein and Fig. 5 in Ryan et al., 2019). Runoff limits and snow lines simulated by MAR are often high, but differences between melt seasons are in qualitative agreement with MODIS observations. On average, $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$, as well as RACMO snow lines, fall below MODIS and variability from year to year appears suppressed.

RACMO's firn module and IMAU-FDM are very similar, apart from the coarser vertical resolution of the former, and for the remainder of the discussion, we focus on IMAU-FDM to establish the main causes for the differences in maxΥ_rcm between MAR and the RACMO family of models.

At the scale of individual melt seasons, daily MAR data shows strong drops in Υ_rcm during cold spells (Fig. 2). IMAU-FDM shows only moderate drops but the smoother curve is due to the coarser 10 d temporal resolution of the data. Sudden drops are not present in MODIS Υ_obs because the actual routing of meltwater is a much slower process than the instantaneous vertical routing in bucket schemes. In slush fields and streams water can flow along the surface for tens of kilometers (Holmes, 1955; Poinar et al., 2015; Yang et al., 2021), at speeds of a few meters per hours in slush (Clerx et al., 2022) or a few kilometers per hour in surface streams (Gleason et al., 2016). Holmes (1955) observed that it took about two weeks after the end of melting before streams ran dry and froze over.

5.2 Why do simulated runoff limits differ?

The substantial differences between runoff limits simulated by MAR and IMAU-FDM (e.g. Figs. 3 and 4) could be caused by (i) differences in RCM simulated accumulation or melt, or (ii) differences in the parameterizations of firn and firn hydrology. A third possible reason are the differences between MAR and IMAU-FDM firn temperatures. We will discuss this aspect in the context of the differences between the firn parameterizations.

On the example of the K-transect we have shown that RCM simulated accumulation and melt (Fig. B1) are generally similar. However, $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ are situated at lower elevations than $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ and because of their lower elevations, melt at $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ is substantially larger than at $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ (Appendix B). Because there can be no runoff above the runoff limit, IMAU-FDM simulated refreezing at $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ is substantially larger than MAR refreezing at $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$. Consequently, we argue that differences in the models' parameterizations of firn and firn hydrology are mainly responsible for the differences between their runoff limits. Thereby the implementation of the bucket scheme plays an important role, namely different choices of (i) irreducible water content, (ii) firn layer depth, (iii) the depth at which runoff can occur, and (iv) the thickness of individual model layers.

IMAU-FDM's large refreezing potential is the main reason for its low runoff limits. The refreezing potential is large due to (i) the relatively low firn temperatures, (ii) the relatively high irreducible water saturation at higher firn densities, and (iii) the thick firn layer (up to 100 m) which offers ample amounts of firn air content in which meltwater can refreeze. IMAU-FDM's condition that runoff can only occur at the bottom of the firn pack, is also responsible for the runoff limit being relatively immobile. Before $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ can propagate to higher elevations, the pore space of the thick firn pack needs to be filled. However, once a grid cell's firn has lost its pore space, this grid cell will nearly always remain runoff area, even during weak melt years: apart from the pore space in the seasonal snow, there is no more possibility to store meltwater. This explains (i) why in IMAU-FDM the uppermost elevation of fully icy firn roughly coincides with $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ (Fig. 7), (ii) why in moderate melt years $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ does not drop substantially below the elevation of fully icy firn, and (iii) why high-elevation melt in extreme melt years cannot run off and instead refreezes, as indicated by the strong firn warming in 2012 (Fig. C3).

RACMO's surface albedo parameterization further contributes to immobilizing $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$. During intense melt seasons, RACMO shows a pronounced step change in surface albedo that coincides with $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ (see Fig. 5c for the situation in the summer of 2012). The higher albedo above that step change reduces melt and also the likelihood of percolation to the bottom of the firn where runoff could take place. Furthermore, reduced melt above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ reduces the amount of water available for refreezing which slows down the loss of firn pore space. The albedo step change is caused by RACMO's ELA coinciding with $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$, a situation which occurred every fourth melt season during the time period 1990–2020. Below the ELA, RACMO albedo is prescribed based on MODIS imagery (see Sect. 3.2.2). Above the ELA, the albedo is calculated based on snow albedo parameterizations independent of MODIS data. MAR does not show discontinuities in albedo, also not in 2012 (Fig. 5c) which is the only melt season where MAR's ELA coincides with $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$. It appears that MAR's albedo parameterization, which does not use remote sensing data, allows for a more smooth transitions of surface albedo across $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$.

MAR's firn temperatures are warmer than IMAU-FDM (Fig. C5), the irreducible water saturation below 1 m depth is smaller than in IMAU-FDM and the simulated firn pack is more shallow reaching only to 21 m depth. This means that MAR's refreezing potential is smaller and allows for stronger fluctuations in $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$, as compared to IMAU-FDM. Runoff in MAR occurs also from areas of porous firn (Fig. 7), which does not occur in IMAU-FDM. The reason is MAR's parameterization which states that $\mathrm{1}/\mathrm{3}$ of meltwater reaching an ice lens runs off immediately while the remaining $\mathrm{2}/\mathrm{3}$ are routed further to depth. This parameterization mimics lateral runoff of meltwater on top of low-permeability ice slabs (MacFerrin et al., 2019) and allows $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ to fluctuate in-between the elevation of depleted firn pore space and the highest elevation where ice layers are simulated in the otherwise porous firn.

It remains unclear why MAR firn temperatures are warmer and show a less smooth spatial distribution than RACMO (e.g. Fig. C3). Spatial discontinuities in MAR firn temperatures were already shown to exist Greenland-wide (Vandecrux et al., 2024). The same publication also shows that MAR firn temperatures are typically higher than for RACMO and hypothesizes this might be linked to numerical instabilities in MAR's firn module. However, MAR's irregular spatial pattern could also be caused by the coarser firn layers and the dynamic vertical discretisation. The latter refers to MAR's merging of adjacent layers of similar properties in order to keep a higher number of layers available to represent the first meter of snow. It can occur that individual MAR pixels have only one layer of ∼20 m in thickness situated below 19 thin layers resolving the first meter of the snowpack. As a result, in some pixels the 10 m depth temperature refers to the temperature of a layer covering a large depth interval, for other pixels to a much thinner layer close to 10 m depth. In IMAU-FDM, the firn is much finer resolved and a comparison to measured firn temperatures at a certain depth (Fig. C5) always compares to a thin model layer very close to that depth. An alternative explanation for the colder IMAU-FDM firn temperatures would be that the Figs. 6, C3 and C5 give a wrong impression because latent heat in IMAU-FDM is released at depths greater than the max. 20 m shown in the figures. If this were the case, then IMAU-FDM depth-integrated firn temperatures would be warmer than the visualized top 20 m. However, this is not the case: During the strongest melt season of 2012, IMAU-FDM meltwater percolation reached a maximum of ∼15 m depth directly above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ and ∼5 m at KAN_U. IMAU-FDM's relatively high irreducible water saturation hinders deep percolation.

5.3 Simulated runoff limits influence total runoff

We find that in intense melt years, MAR simulates up to 29 % more runoff than IMAU-FDM along the K-Transect. This difference is mainly due to MAR runoff from above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$. All $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ are located further inland than $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ and in the year 2012, the distance between the two runoff limits reaches ∼75 km (Fig. 5a). While average MAR runoff between the two runoff limits is modest compared to average runoff over the RCM's common runoff area (Fig. 5e), the considerable distance causes total MAR runoff from above $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ to become relatively large. In melt seasons of intermediate intensity, MAR and IMAU-FDM maxΥ_rcm are located closer to each other (Fig. 5b) and total runoff between them is relatively small.

Although 2012 and 2019 appear as outliers when compared to most other melt seasons, the trend towards larger differences in strong melt seasons is a consequence of $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ varying weakly with melt intensity while $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ fluctuates strongly. The ice sheet hypsometry amplifies this effect. As the ice sheet surface becomes increasingly flatter towards higher elevations (van As et al., 2017), elevation differences between the two maxΥ_rcm translate into large horizontal offsets. In strong melt seasons, $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ are located at elevations where the surface slope is shallow and horizontal distance to the lower $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ becomes large (e.g. Fig. 5a and b).

5.4 Implications

Our analysis focuses on the K-transect, which is located where differences in maxΥ_rcm are at their maximum (Fig. 3). However, other studies indicate our findings are valid elsewhere on the ice sheet. Spatial discontinuities in MAR firn temperatures were already shown to exist Greenland-wide by Vandecrux et al. (2024). Tedstone and Machguth (2022) focused on firn areas that experience surface runoff and found that 1985–2020 MAR and RACMO simulated cumulative runoff above a certain reference elevation differ by a factor of two. Given the relationship shown in Fig. 8 and our explanation why the difference between the two maxΥ_rcm increases with melt season intensity, we expect runoff limits to diverge further in a warmer future climate. Indeed, Glaude et al. (2024) show that by the year 2100, under identical SSP5-8.5 high emissions forcing, the runoff limits of RACMO and MAR differ strongly over most of the ice sheet. The consequence is a twofold larger simulated annual surface mass loss in MAR than in RACMO (Glaude et al., 2024).

Uncertainty in future Greenland surface mass balance will grow with continued warming, and uncertainties in simulating Greenland's firn area contribute strongly to overall uncertainty. As both models demonstrate strengths and weaknesses in reproducing MODIS Υ_obs and maxΥ_obs, it is unknown which simulates total runoff more accurately. Nevertheless, combining the strengths of the models might be a first step to improve the simulation of the surface mass balance of Greenland's firn area.

RACMO, and consequently also IMAU-FDM, might benefit from a revised bare-ice albedo parameterization. The existing parameterization leads to step-like changes in albedo at the runoff limit during intense melt seasons. IMAU-FDM simulates a finely resolved and deep firn column, but this leads to a relatively immobile runoff limit when combined with a standard bucket scheme where runoff can take place only at the base of the firn. In a first step, IMAU-FDM could include a parameterization that mimics lateral runoff whenever percolating water encounters an ice layer, akin to the parameterization included in MAR.

Once the potential numerical instabilities in MAR's firn simulation are resolved, the model might benefit from simulating a deeper and more finely resolved firn column. The current coarse resolution and the merging of layers impede comparisons to measurements and challenge assessment of model performance. MAR includes a parameterization mimicking the effect of ice slabs on runoff. However, we assume that the large differences between $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ and maxΥ_obs are related to this parameterization and we suggest that it should be calibrated. The minimum thickness, required for an ice layer to trigger runoff, could be set based on Jullien et al. (2025) who provide first empirical evidence for the minimum ice slab thickness supporting lateral runoff. Altering the runoff ratio from $\mathrm{1}/\mathrm{3}$ to another value would not directly influence $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$, but controls how much water percolates to depth and thus influences refreezing and firn structure, such as the formation or thickening of ice layers.

Beyond these initial modifications, the models could replace the bucket scheme with more physical simulations of snow and firn as applied by Wever et al. (2014, 2016); Langen et al. (2017); Vandecrux et al. (2020). Besides the inclusion of preferential percolation, these approaches also allow for temporary storage of meltwater in snow and firn, which plays an important role in shaping firn structure. Observations since 2012 at the KAN_U site show that the ice slab is not getting buried as simulated by IMAU-FDM. Instead, the depth of the slab remained roughly constant (Fig. C4a). The ice slabs are of low permeability which causes meltwater to pond in slush at their surface (Clerx et al., 2022) and to refreeze partially, over the course of a melt season, as superimposed ice (Tedstone et al., 2025). This mechanism, by which ice slabs mainly thicken, is absent in an instantaneous bucket scheme. Both RCMs currently do not permit slush formation and even thick ice layers must remain “permeable” for meltwater to be routed vertically. Removing these constraints by adopting more physical firn simulations might improve the models' representation of melting firn.

6 Conclusions

We developed a flexible method to detect visible runoff limits from MODIS and compared the results to modelled runoff limits from IMAU-FDM and MAR. We found large differences not only between remotely sensed and modelled data, but also between the two models. IMAU-FDM simulated runoff limits are on average somewhat lower than MODIS, and variability from year to year is strongly reduced. On average, MAR simulates substantially higher runoff limits than MODIS, but the magnitude of yearly fluctuations of MAR's runoff limits are mostly similar to MODIS. Both MAR and IMAU-FDM use a bucket scheme that routes water vertically through the firn. Differences in the implementation of the bucket schemes are an important reason for the deviations between MAR and IMAU-FDM runoff limits: (i) in MAR a fraction of the meltwater runs off when it encounters an ice layer inside the firn, (ii) the amount of pore space and cold content varies between the two models because they simulate different firn depths, and (iii) IMAU-FDM allows for a higher irreducible water saturation. Furthermore, the firn layer in MAR is generally warmer which reduces the retention capacity and promotes runoff. It is assumed that MAR's warmer firn layer is caused by currently unresolved instabilities in the models firn module.

We compare total simulated RCM runoff along the K-transect and we find that MAR total runoff exceeds IMAU-FDM by up to 29 %. We show that in strong melt seasons MAR and IMAU-FDM runoff limits are separated by large horizontal distances, which is the main reason for the difference in total runoff. Any differences in RCM ablation area runoff are eclipsed by the amount of runoff that MAR simulates, in strong melt years, above the IMAU-FDM runoff limit. Ice sheet hypsometry contributes to the large horizontal distance between the two runoff limits: the ice sheet surface slope becomes increasingly shallow with altitude and relatively small differences in runoff limit elevations translate into large horizontal distances.

Increased melting is anticipated for the future. This means the situation where the two models diverge the most will become more frequent. We hypothesize that simulated runoff will further diverge and uncertainty will grow. We conclude that newly formed runoff areas will play a major role in Greenland's future mass balance. Reliably simulating the surface mass balance of melting firn is key to faithfully anticipate the future of the Greenland Ice Sheet.

Appendix A: Calculation of Υ_obs along flowlines

A1 Calculation along flowlines

We create polygons by calculating flowlines which are buffered by $p_{\mathrm{w}}/\mathrm{2}$. This approach creates polygons of arbitrary shape and direction (Fig. A1), here termed flowline-polygons. Even in complex topography, the direction of the flowline-polygons is always roughly perpendicular to the surface slope.

Figure A1 Flowlines (orange) and flowline polygons (blue shaded areas) at Sermeq Kujalleq (Jakobshavn Isbræ). Darker shades of blue indicate overlapping polygons.

We chose to calculate flowlines based on surface velocity fields rather than surface slope (cf. Machguth and Huss, 2014). The advantage is a straightforward algorithm, as described in the following. We calculate flowlines following Fig. 3 in Cabral and Leedom (1993), using Greenland ice sheet surface velocity fields in x and y direction. Our algorithm starts at seed-points and then progresses downhill from gridcell to gridcell. A flowline enters a cell at a certain point along its margins and based on entry point, flow direction within the cell and cell size, the algorithm then calculates the point where the flowline leaves the cell and enters the following cell. A flowline ends when it reaches the ice sheet margin.

There are cases where flow directions of neighboring cells are conflicting and the algorithm would send the flowline immediately back to the cell where it came from. Such conflicts are solved by calculating the average flow direction of the two grid cells in question. The flow line then continues in average flow direction through one of the two cells.

Seed-points are created by first drawing a polygon that follows roughly the 2400 m a.s.l. elevation contour in the south of the ice sheet and descends towards the 1800 m a.s.l. contour in the north. Along the polygon, seed-points are created automatically every 15 km. Eventually, all flowlines are buffered by $p_{\mathrm{w}}/\mathrm{2}=\mathrm{10}$ km to create flowline-polygons. Given the width of the flowline-polygons (p_w=20 km) and 15 km spacing of the seed points, a certain overlap of the polygons occurs and is wanted (Fig. A1). More closely spaced polygons provide a higher spatial resolution of Υ_obs and make it easier to detect outliers. On outlet glaciers polygons overlap due to confluence (Fig. A1). There are also cases where polygons overlap for most of their length due to a combination of specific flow patterns and location of the seed points. The polygons were sifted manually to remove such polygons. The result is a set of 510 flowline-polygons (see Fig. 1).

A2 Accounting for background spatial variability of albedo

Our algorithm uses daily MODIS MOD10A1 albedo maps to assess spatial variability of albedo σ_α. MODIS records changes in α and σ_α as surface characteristics and hydrology evolve over the duration of a melt season. However, the satellite images also capture pattern in α that are persistent in space and time. Such persistent albedo features typically originate from topographic undulations or rock outcrops. Where persistent albedo features are frequent, they impact σ_α and interfere with detecting Υ_obs. The original approach by Machguth et al. (2023) did not include any correction for the potential impact of persistent albedo features on Υ_obs. The updated approach used here now includes a correction as described in the following.

We calculate a Greenland-wide map of background σ_α, based on daily arrays of σ_α from before the start of the melt season. (i) From each spring of the 22 years 2000 to 2021, 20 daily arrays of σ_α are selected. (ii) We then calculate grid cell values of an initial background σ_α array as the median of up to 440 (22 years × 20 d) daily values (the actual number of data points is smaller due to frequent clouds or data issues). The large north-south extent of Greenland requires to vary the 20 d time-window across latitudes. Up to ∼75.5° N the time window are the days of year (DoY) 110–130, between ∼75.5 and ∼80° N DoY 120–140, and north of ∼80° N DoY 130–150. (iii) The final array of background σ_α is calculated by subtracting the mean of all grid cells, calculated from the initial background σ_α array, from each grid cell. Any resulting negative values are replaced by zero.

In detecting daily Υ_obs, the final array of background σ_α is subtracted from every daily array of σ_α. The thresholds for σ_α, used in the original algorithm by Machguth et al. (2023), remain unchanged as the background σ_α array consists mostly (82 %) of zeros.

A3 Modified filtering for outliers

Candidates for Υ_obs require filtering to remove false positives (Machguth et al., 2023). We apply the same automated approach in two stages but the filtering of the last valid candidates has been simplified (Sect. 4.4 in Machguth et al., 2023). If a suspicious last candidate is detected, then the updated algorithm searches for valid detections within a time window of ±6 d and a circle of 75 km. The suspicious candidate is labeled invalid if it exceeds the median elevation of all nearby valid detections by >75 m. If the number of nearby valid detections is too small to calculate a median, the suspicious candidate is labeled “valid”. The number of removed candidates remains similar under the updated filter algorithm, but there is no more risk of consulting distant Υ_obs when evaluating reliability of candidates.

A4 Comparison to Landsat-derived visible runoff limits

We compared MODIS Υ_obs to annual maxima of Landsat visible runoff limits RL, using annual maximum RL at 1 km posting (see methods in Tedstone and Machguth, 2022). We first iterated through each flowline polygon, identifying all the Landsat RL which fall inside it, then generated median Landsat RL for all data in that polygon on a particular day. We only compare MODIS and Landsat on days when retrievals were made by both approaches and comparisons were only done for those flowline polygons located in areas for which Tedstone and Machguth (2022) applied their Landsat algorithm. Among smaller excluded areas on the west coast and in the north, no comparison was possible for the entire east coast south of ∼76° N.

The comparison is shown in Fig. A2 and yields a linear regression that falls very close to the line of identity. The bias between the two datasets is small, on average MODIS Υ_obs falls 26 m below Landsat RL. The comparison yields R²=0.81, which is somewhat lower than the R²=0.87 of the evaluation of the Machguth et al. (2023) algorithm against Landsat visible runoff limits. However, the comparison in Machguth et al. (2023) was restricted to the west coast which is the area where MODIS and Landsat visible runoff limits are most reliable. Furthermore, the comparison shown in Fig. A2 focuses on Landsat annual maximum RL while Machguth et al. (2023) used all individual Landsat visible runoff limit retrievals followed by detection and removal of likely erroneous Landsat visible runoff limits. Here we do not apply any cleaning to the Landsat RL.

Figure A2 Comparison of daily MODIS Υ_obs and Landsat derived visible runoff limits (RL; Tedstone and Machguth, 2022). The number of samples is n=3880.

Qualitatively, we conclude that the improved MODIS algorithm compares similarly to Landsat RL as did the original MODIS algorithm by Machguth et al. (2023). The latter, however, was restricted in its applicability to the western flank of the ice sheet. We find the largest deviations between the improved MODIS algorithm and Landsat at the north-eastern flank of the ice sheet. For example, the point cloud located below the line of identity at Υ_obs≈850 m a.s.l. (see Fig. A2) concerns MODIS and Landsat retrievals from the vicinity of flowline NE (Fig. 1).

Appendix B: Differences in accumulation and melt close to $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ and $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$

We explore differences in melt M_rcm and accumulation C_rcm at maxΥ_rcm and investigate whether they could explain the differences in modelled runoff limits. For clarity, we focus the analysis on the K-Transect. First we compare annual accumulation sums in RACMO (C_RACMO) and MAR (C_MAR). We sum up C_rcm over hydrological years (1 September to 31 August) and average over a zone that encompasses all annual maxΥ_rcm of IMAU-FDM and MAR. We focus on this zone rather than the entire K-Transect as we want to examine differences close to the maxΥ_rcm. We observe a high correlation of annual accumulation simulated by the two RCMs ($C_{\text{RACMO}}=\mathrm{0.09}+\mathrm{0.93}{C}_{\text{MAR}}$; R²=0.92, p<0.001). Average C_RACMO (0.44 ± 0.08 m w.e.) exceeds average C_MAR (0.37 ± 0.09 m w.e.). Next we regress annual $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ vs. C_MAR and $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ vs. C_RACMO. Both regressions do not yield statistically significant relationships, indicating that differences in C_rcm cannot explain the differences between the models' runoff limits.

Figure B1 Comparison of seasonal simulated melt in MAR and RACMO along the K-transect. Melt is averaged over the grid cells located in-between the lowest and highest of all Υ_rcm and over 1 June to 30 September of each year. (a) Linear regression of MAR and RACMO seasonal melt. (b) Scatterplot of seasonal melt in MAR and RACMO vs. ${\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ and ${\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$, respectively. In both subplots the 2012 and 2017 melt seasons are marked with a plus sign.

Second, we compare melt for the same zone and summed up over each melt season, defined as 1 June to 31 August. We find that M_RACMO and M_MAR are highly correlated but RACMO melt is biased low in comparison to MAR (Fig. B1a). However, the bias is small or close to zero for moderate and low melt seasons, respectively. The differences in M_rcm might be explained by RACMO having on average a higher surface albedo (0.79 ± 0.02) as MAR (0.77 ± 0.02). Regressing annual maxΥ_rcm against M_rcm reveals a stark contrast between the two RCMs (Fig. B1b). For a given amount of melt, $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ exceeds $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ by up to ∼450 m. The $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ show a weak dependency on M_RACMO while $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ depend more strongly on M_MAR. Differences between M_RACMO and M_MAR apparently cannot explain the large differences in Υ_rcm either.

Third, we compare C_rcm and M_rcm simulated at the RCM grid cells that coincide with each annual Υ_rcm. We find rather similar average C_RACMO at $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ (0.40 ± 0.07 m w.e.) and C_MAR at $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ (0.37 ± 0.09 m w.e.). Average M_rcm at maxΥ_rcm is higher in RACMO (0.59 ± 0.21 m w.e.) than in MAR (0.34 ± 0.12 m w.e.). This is consistent with the above established low bias of M_RACMO because the $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{IMAU-FDM}}$ are located at substantially lower elevations where melt is higher. The comparison of C_rcm and M_rcm at annual maxΥ_rcm reveals an important difference between the models: in IMAU-FDM, the runoff limit is typically located where summer melt exceeds annual accumulation ($C_{\text{RACMO}}-M_{\text{RACMO}}=-\mathrm{0.19}$ ± 0.25 m w.e.); in MAR melt and accumulation at $\text{max}{\mathrm{{\rm Y}}}_{\text{rcm}}^{\text{MAR}}$ are similar ($C_{\text{MAR}}-M_{\text{MAR}}=\mathrm{0.03}$ ± 0.14 m w.e.).

Appendix C: Additional Figures

Figure C1 Seasonal evolution of Υ_rcm simulated by MAR and IMAU-FDM, as well as Υ_obs detected from MODIS. The comparison is shown for a flowline-polygon located at around 66° N on the west coast (region SW, see Fig. 1).

Figure C2 Seasonal evolution of Υ_rcm simulated by MAR and IMAU-FDM, as well as Υ_obs detected from MODIS. The comparison is shown for a flowline-polygon located at around 80° N (region N, see Fig. 1).

Figure C3 Comparison of RCM simulated 10 m firn temperatures along the K-transect, 2010 to 2020. Data to the left are simulated by IMAU-FDM; MAR data are shown to the right.

Figure C4 Evolution 1980–2020 of RCM simulated firn density ρ in the vicinity of the KAN_U site (K-transect at 1840 m a.s.l.). Dotted areas show where ρ>830 kg m⁻³, i.e. exceeding pore close-off density. Green dots mark in situ measured depths of the top of the ice slab.

Figure C5 Comparison of RCM simulated and measured firn temperatures at the KAN_U site (1840 m a.s.l.) and for the years 1980 to 2020. (a) Firn temperatures for the top 20 m simulated by FDM; (b) top 20 m firn temperatures modelled by MAR; (c) comparison of modelled and measured firn temperatures at 10 m depth (Charalampidis et al., 2016; How et al., 2022; Vandecrux et al., 2024; Vandecrux, 2023). White dots in subplots (a) and (b) denote the top of the ice slab surface according to the measurements summarized in Rennermalm et al. (2021).