Assessing Bare Ice Albedo Simulated by MAR over the Greenland Ice Sheet (2000-2021) and Implications for Meltwater Production Estimates

. Surface mass loss from the Greenland ice sheet (GrIS) has accelerated over the past 10 decades, mainly due to enhanced surface melting and liquid water runoff in response to atmospheric warming. A large portion of runoff from the GrIS originates from exposure of the darker bare ice in the ablation zone when the overlying snow melts, where surface albedo plays a critical role in modulating the energy available for melting. In this regard, it is imperative to understand the processes governing albedo variability to accurately project future mass loss from 15 the GrIS. Bare ice albedo is spatially and temporally variable and contingent on non-linear feedbacks and the presence of light-absorbing constituents. An assessment of models aiming at simulating albedo variability and associated impacts on meltwater production is crucial for improving our understanding of the processes governing these feedbacks and, in turn, surface mass loss from Greenland. Here, we report the results of a comparison of the bare ice extent and albedo 20 simulated by the regional climate model Modèle Atmosphérique Régional (MAR) with satellite imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS) for the GrIS below 70° N. Our findings suggest that MAR overestimates bare ice albedo by 22.8% on average in this area during the 2000-2021 period with respect to the estimates obtained from MODIS. Using an energy balance model to parameterize meltwater production, we find this bare ice albedo bias can 25 lead to an underestimation of total meltwater production from the bare ice zone below 70° N of 42.8% during the summers of 2000-2021.


Introduction
Global mean sea level (GMSL) rise has significantly accelerated over the past decades (Chen et al., 2017), in part as a consequence of the acceleration in Greenland ice mass loss (Aschwanden et al., 2019). Ice mass loss from the Greenland ice sheet (GrIS) was one of the largest contributors to GMSL rise in the period 1901-2018 with 17-32% and will likely remain so by the 5 end of this century (Fox-Kemper et al., in press). According to Fox-Kemper et al. (in press), the GrIS' contribution to GMSL will constitute ~17% by 2100 for the Shared Socioeconomic Pathway SSP5-8.5 (Riahi et al., 2017). From this point of view, it is imperative to improve model representation of physical processes responsible for ice mass loss to better constrain projections of the future contribution of the GrIS. In this regard, evaluation of climate model outputs vs. 10 observational data can provide insight into the model's ability to represent the physical processes at play and can subsequently highlight regions for model improvement .
The total mass loss from Greenland can be separated into surface losses (e.g., runoff) and frontal losses at the terminus of outlet glaciers (e.g., calving). For the 2000-2018 period, 55% of 15 Greenland's mass loss originated from surface mass balance (SMB; the balance between accumulation and ablation at the ice sheet surface) and 45% from ice discharge from outlet glaciers (Mouginot et al., 2019). The SMB losses from the GrIS have been increasing since the late 1990s, driven primarily by an increase in melt and subsequent liquid water runoff in response to recent atmospheric warming . An increase in summer surface air 20 temperatures of ~+2° C has also been observed over Greenland since the early 1990s (Hanna et al., 2012;Box, 2013). This has increased runoff by 40%, while the contribution of changes in precipitation, sublimation, and erosion since the early 1990s are not as substantial (van den Broeke et al., 2016. The recent increase in observed runoff has also been suggested to be linked to changes in atmospheric circulations around Greenland (Hanna et al., 2014;Tedesco and Fettweis, 25 2020). (Hanna et al., 2018;McLeod and Mote, 2016) suggest that a significant increase in high summer pressure blocking over Greenland (Greenland Blocking Index) since the 1990s has been a major driver of the recent increase in surface melt over the GrIS (Fettweis et al., 2013;Tedesco et al., 2016a). Changes in atmospheric circulation promoting enhanced runoff have been characterized by increased shortwave incoming solar radiation . This, in turn, can increase the absorbed solar radiation, depending on surface albedo, which can be summarized as the ratio of the energy reflected off a unit area of material over the energy incident over that area.
Besides temperature, albedo also strongly controls surface melting and runoff over the GrIS. More specifically, broadband albedo refers to the wavelength-dependent albedo integrated 5 over the full spectral range, weighted by the contribution of each wavelength. In the case of solar radiation, most of the contribution to broadband albedo comes from the visible wavelengths (300-700nm), since most of the solar energy is concentrated within this band (Liang and Wang, 2020).
Albedo is one of the key players in the energy balance for the bare ice of the GrIS, which is exposed when the overlying seasonal snow melts. Snow is characterized by a high albedo of ~0.7-0.85 10 (Alexander et al., 2014), while bare ice is compacted, densified and aged snow (Wiscombe and Warren, 1980) and is characterized by a low albedo of 0.31-0.57 (Wehrlé et al., 2021). Bare ice thus absorbs more solar radiation than snow, increasing the energy available for melting. Even though the bare ice zone encompasses only a small fraction (12 ± 2 %) of the GrIS in summer along the margins of the ice sheet (Ryan et al., 2019), the bare ice zone was responsible for 78% 15 of the runoff from the GrIS in the period 1960-2014 (Steger et al., 2017). Since bare ice albedo strongly controls the amount of runoff from the bare ice zone, it is of key importance in controlling GrIS-wide runoff (Tedesco et al., 2008;van Angelen et al., 2012;Alexander et al., 2019).
Bare ice albedo is spatio-temporally variable at different scales in response to non-linear positive feedbacks between absorbed shortwave radiation and surface melt Ryan 20 et al., 2019). Therefore, meltwater production does not only depend on timing and persistence of bare ice exposure but also on other modulating factors, such as snowfall, which can cover bare ice with a bright, highly reflective layer of fresh snow. The appearance of bare ice reduces the overall GrIS albedo, leading to more melting from the bare ice zone as well as feeding a positive feedback mechanism which ultimately leads to an acceleration of surface melting (Tedesco et al., 2011). 25 Bare ice exposure is often associated with the presence of light-absorbing constituents (LACs) on the ice, such as dust, black carbon and organic material (Tedstone et al., 2017). These LACs can reduce surface albedo to as low as ~0.1, (Wientjes et al., 2012;Tedstone et al., 2017Tedstone et al., , 2020, and can subsequently increase melting. Dark bands appear in the bare ice as a consequence of outcropping of ice layers mixed with dust that were deposited in the accumulation zone during the late Pleistocene and Holocene and, later, transported to the lower ablation zone through ice flow (MacGregor et al., 2020). Black carbon and cryoconite, small cylindrical holes (of a few centimeters to a few meters) in the ice surface containing impurities, have also been found to reduce the albedo of bare ice (Cook et al., 2016;Goelles and Bøggild, 2017). (Wang et al., 2018;5 Stibal et al., 2017) found an abundant presence of supraglacial ice algal blooms in the bare ice zone in the Southwestern GrIS, with a direct link between mineral phosphorus in the ice surface and glacier ice algae biomass (McCutcheon et al., 2021). These light-absorbing constituents reduce the bare ice albedo, further enhancing meltwater production and runoff (Tedesco et al., 2016b;Williamson et al., 2018Williamson et al., , 2020Cook et al., 2020). The difficulty in representing bare ice albedo in 10 climate models partly originates from a lack of understanding of LACs and may result in a reduced accuracy of runoff projections (Alexander et al., 2014).
In this study, we evaluate the performance of the Modèle Atmosphérique Régional (MAR), a regional climate model especially developed for simulating polar climates , in simulating bare ice extent, bare ice albedo and meltwater production by comparing MAR's 15 model output with satellite imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS). This study complements the study by (Alexander et al., 2014) who focused on the GrISwide albedo. Here, we specifically focus on the bare ice zone below 70° N, which is currently responsible for the majority of meltwater production from the GrIS (Steger et al., 2017). We evaluate MAR on a range of spatial resolutions during June, July and August in 2000-2021. We 20 use an energy balance model to parameterize the meltwater production and to analyze the effect of a bias in observed and modeled bare ice albedo on estimates of meltwater production.

The MAR RCM
In this study, we use the Modèle Atmosphérique Régional (MAR) version 3.12 regional 25 climate model (RCM), which simulates the coupled surface-atmosphere system over the Greenland region (Gallée, 1997;Ridder and Schayes, 1997;Lefebre, 2003;Fettweis et al., 2017) and is forced by reanalysis data or climate model output. Specifically, we force MAR at the lateral boundaries and ocean surface with 6-hourly reanalysis output from ERA5 (Hersbach et al., 2020), produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). The MARv3.5.2 is discussed and validated over the GrIS , with updates to MARv3.11 discussed in (Fettweis et al., 2021). We point out that in the version we use in this study, the geographical projection has been changed to the Standard Polar Stereographic EPSG 3413 from a previously used custom projection. An issue within the code impacting the snow 5 temperature at the base of the snowpack has also been corrected in MARv3.12 and it now includes a continuous conversion from rainfall to snowfall from 0° C to -2° C as input for the snow model instead of a fixed value of -1° C (Fettweis, personal comm.). The atmospheric component of the model is described by (Gallée and Schayes, 1994) and the Soil Ice Snow Vegetation Atmosphere Transfer scheme (SISVAT) is used as the surface component of the model (Ridder and Schayes, 10 1997). The surface model incorporates the snow model CROCUS (Brun et al., 1992), which simulates a set number of layers of snow, ice, or firn with variable thickness and transports energy and mass between each layer. The CROCUS model also provides snow grain properties, which are used to simulate snow albedo. In this study, we run the MAR model over Greenland and produce daily output of variables pertaining to the atmosphere and ice sheet surface in this region 15 at horizontal spatial resolutions of 6.5, 10, 15, and 20 km.

MODIS data
We obtained MOD09GA Version 6 (Vermote and Wolfe 2015) surface reflectance images over the GrIS from the Moderate Resolution Imaging Spectroradiometer (MODIS) on board NASA's Terra satellite through Google Earth Engine (Gorelick et al. 2017). We use daily summer 20 (June, July, and August; JJA) images over the period 2000-2021 with a horizontal spatial resolution of 500 m. Corrections have been applied to this product for atmospheric conditions such as aerosols, gasses and Rayleigh scattering (Vermote and Wolfe 2015). We also collected daily snow cover images from MOD10A1 Version 6 from Google Earth Engine over the same period.
This product includes a daily cloud mask which we use in this study to flag clouds in the 25 MOD09GA images. The MOD10A1 product also contains daily broadband albedo values (Hall et al. 2016), though albedo values above a latitude of 70° N may be positively biased due to viewing geometry and large solar zenith angle (Alexander et al., 2014). Consequently, we omit MOD10A1 albedo data above 70º N from our albedo and meltwater production analysis. Note that we do include MOD09GA surface reflectance data above 70º N in our bare ice extent analysis. We choose the summer 2000-2021 study period to accommodate the observation period of MODIS and to account for the seasonal variability of bare ice exposure on the GrIS, when surface albedo has the largest impact on SMB (Alexander et al., 2014). To allow for a daily pixel-by-pixel comparison, we first use GDAL (GDAL/OGR contributors) to reproject the daily MODIS data (MOD09GA and MOD10A1) to the MAR's native projection and simultaneously rescale it to the resolution of 5 each of the MAR products.

Bare ice extent
Bare ice is exposed when the snow blanketing it is removed through surface melting. Most of the areas where bare ice occurs are located in the ablation zone, where ablation is larger than accumulation and the SMB is negative. At the transition between the ablation and accumulation 10 zones lies the equilibrium line altitude (ELA), denoting the elevation where ablation is equal to accumulation and the SMB is 0 . In order to study the behavior of bare ice exposure we determine a long-term average ELA from the daily MAR outputs of SMB over the GrIS. We estimate the average ELA at 1679 m a.s.l. for the period 2000-2021 over the entire ice sheet as the 95th percentile value of the elevation values in the ablation zone. Taking the 95th 15 percentile of the long-term average values supports the omission of sporadically high ablation cell detections and provides a conservative estimate of the ELA. We note that this method may provide a conservative estimate of the bare ice extent during warm, high-melt years, as the ELA in such years may lie at a higher elevation than the long-term average ELA. Then, we constrain bare ice as modeled by MAR to cells below the long-term average ELA. This is a first simplified estimate 20 of the bare ice extent, which is further refined by the following two conditions: 1) snow is absent (i.e. snow depth is 0 m) and 2) the average density in the top 1 m exceeds 907 kg/m 3 . A thin layer of fresh snow (300 kg/m 3 in MAR) could cover the ice (920 kg/m 3 in MAR) following a brief snowfall event. Solar radiation will not attenuate much in a thin layer of snow and will thus not significantly affect absorption into the underlying ice (Warren et al., 2006). A thin layer of fresh 25 snow will lower the density of the top layer, however. Therefore, setting a lower limit of 907 kg/m 3 for the average density of the top 1 m allows for 2 cm of fresh snow to cover the ice, while also allowing the cell to still be detected as bare ice and not as snow. Taking the average density also ensures that ice lenses are not detected as bare ice. We use the static ice mask and digital elevation model (DEM) of the GrIS as described by the Greenland Ice Mapping Project (GIMP) to select areas where the ice sheet is present (vs. where land is present) and to produce the satellite-derived extent and elevation of the GrIS (Howat et al., 2014).
We note that the static GIMP ice mask and DEM are constructed from Landsat-7 and RADARSAT-1 imagery acquired between 1999 and 2002, which only overlaps for 3 out of the 22 years of the study period in this study. However, the impact on the estimated bare ice extent should 5 be small given estimates for ice margin retreat rates and thinning rates (Helm et al., 2014;Lesnek et al., 2020;Young et al., 2021). We therefore believe that the static GIMP ice mask suffices for the purposes of this study.
We extract the satellite-derived bare ice extent (BIE) on the GrIS by applying an upper threshold of 0.6 to band 2 (841-876 nm) in the MOD09GA product (Shimada et al., 2016). Following the 10 same study, we define pixels with reflectance values above 0.6 in band 2 as snow. We define the annual maximum BIE as the area covered by those pixels that are detected as bare ice for a minimum of 10% of the observed days in JJA in one year. The aim of this is providing a conservative estimate of bare ice extent while ensuring omission of sporadic and erroneous bare ice detections by MODIS, such as superimposed ice and meltwater lakes and streams. We define 15 a lower estimate of the MODIS-derived annual maximum BIE as the area covered by the pixels below the long-term average ELA that are detected as bare ice for a minimum of 10% of the observed days in JJA. We also define an upper estimate of the annual maximum BIE, which includes: 1) the area found for lower estimate and 2) the area covered by the pixels that are flagged as clouds in MOD10A1 for a minimum of 90% of the observed days in JJA. For pixels that are 20 covered by clouds for more than 90% of the observed days in JJA, the view of the surface of the GrIS is obstructed to such an extended degree that bare ice cannot be detected for more than 10% of the observed days in JJA. This automatically excludes them as a possible bare ice pixel, leading to potentially missed bare ice area.
We use forecast verification to quantify MAR's ability to simulate bare ice vs. snow. We Index. This statistic indicates a perfect forecaster with a score of 1, an underforecaster with a score lower than 1, and an overforecaster with a score higher than 1 (Wilks, 2011).

Bare ice albedo
We evaluate the performance of MAR in simulating bare ice albedo by comparing MAR's modeled albedo values with the albedo values observed by MODIS on the overlapping BIE. The 5 output from MAR contains daily albedo values over the entire GrIS. The bare ice albedo scheme in MARv2 originally consisted of simply assigning a fixed value of 0.55 to bare ice albedo.
Nevertheless, the improved MARv3 we use in this study simulates bare ice albedo as a function of accumulated surface water height and slope of the ice sheet, following an exponential relation between pure bare ice albedo and water albedo (Alexander et al., 2014). MAR includes lower and 10 upper boundaries for the bare ice albedo of 0.5 and 0.55. However, every year, large areas with albedo values lower than 0.5 are observed by MODIS and by some PROMICE automatic weather stations (Tedesco et al., 2016b). An analysis by (Wehrlé et al., 2021) shows that the average albedo from the 20 stations included in the study is lower than 0.5 for more than a month during the melt season. Therefore, values for the surface albedo below 0.5 can be considered a common event 15 across the bare ice zone. Such low albedo values in part result from the presence of LACs on the bare ice, which are not taken into consideration in the MARv3 bare ice albedo scheme. Low albedo values can also arise from accumulated meltwater on the surface of the ice in the form of streams and lakes, but the relative effect of meltwater has been estimated to be smaller than that of impurities (Ryan et al., 2018).

Meltwater production
We use an energy balance model to parameterize the energy available for meltwater production over the bare ice zone and to isolate the effect of albedo on meltwater production estimates from the bare ice zone below 70° N. Following (Pellicciotti et al., 2008), we parameterize the energy available for meltwater production as: with daily values for meltwater production (ME), albedo ( ), downward shortwave radiation (SWdown), net longwave radiation (LWnet), and sensible and latent heat fluxes (SHF and LHF) simulated by MAR. The numerator on the right hand side is equal to the energy available for melt.
Dividing by density ( w = 1000 kg/m 3 ) and the latent heat of fusion of water (Lm = 3.34⋅10 5 J/kg) gives the potential meltwater production in mmWE/day. For this purpose, we use the MAR outputs generated at a horizontal spatial resolution of 6.5 km. We determine the parameters a, b, c, and d by finding the minimum of this unconstrained multivariable function on a daily basis. Figure 1   5 shows the linear regression between the parameterized meltwater production and the meltwater production simulated by MAR, with an R 2 of 0.92, mean bias of -0.728 mmWE/day and root mean square error of 3.97 mmWE/day. As seen in Figure 1 and from the negative mean bias, the parameterization tends to underestimate meltwater production slightly relative to the meltwater simulated by MAR. This could be due to the fact that MAR calculates meltwater production every 10 minute and the parameterization calculates melt only once per day since it uses daily MAR output.
Moreover, the feedbacks between air and surface processes are not captured in the parametrization scheme. Lastly, days with melt occurring only during a part of the day occur at the beginning and end of the melt season with the parameters (a, b, c, and d) not being able to fully account for this variation. The fractional contribution to meltwater production of each constituent in Equation 1 15 are calculated by multiplying each parameter with the respective net energy flux and dividing by the total meltwater production on a daily basis.

Figure 1: Linear regression of daily meltwater production simulated by MAR and
parameterized meltwater production using modeled albedo.
We calculate daily meltwater production estimates with the meltwater production parameterization twice, using the same values for the coefficients, once using the albedo modeled by MAR and once using the albedo observed by MODIS. The goal of this is to isolate the effect 5 of bare ice albedo on meltwater production. As a reminder, we exclude cells above a latitude of 70° N to account for the potentially reduced accuracy of albedo values in the MOD10A1 product in this region. In order to increase the fairness of the comparison we include only those areas and days where we simultaneously detect bare ice with both MAR and MODIS. Since we are interested in the effect of bare ice albedo on meltwater production, we only include cells that are melting as 10 prescribed by MAR (>1 mmWE/day). This ensures that any absorbed energy fluxes predominantly go into the enthalpy of fusion of ice, i.e. leading to melt, and not into changing the temperature of the ice. 15 The average number of days when bare ice is exposed during June, July, and August (JJA) for our study period (2000-2021) obtained from MODIS (Fig. 2a) and MAR (Fig. 2b) as well as their difference (Fig. 2c) are shown in Figure 2. The number of bare ice days observed by MODIS

Bare ice extent
show slightly more inland variation, especially in the northern regions of the ice sheet. For instance, a large round feature in the northeast reveals that bare ice is exposed for up to 12 days on 20 average during JJA. The geothermal heat flux map of Greenland created by (Martos et al., 2018) shows a round feature of similar size in the same area with an enhanced heat flux. Increased heat flux from the bedrock to the ice sheet surface could lead to higher ice sheet surface temperatures, which enhances bare ice exposure. This pattern is not captured by the MAR simulations.  Years with a high maximum BIE, such as 2012 and 2019, correspond to known high-melt years . This is as expected since warmer temperatures or positive energy balance anomalies lead to more snowmelt, exposing the underlying bare ice. Additionally, exposed bare ice leads to increased absorption of solar radiation, generating higher melt rates (Ryan et al., 2019). In these high-BIE years MAR overestimates the BIE relative to the 5 observations. This indicates that MAR potentially overestimates the amount of snow that is melted away in these years and exposes more bare ice than is actually the case. The observed and modeled seasonal BIE exhibit a peak from mid-July through mid-August ( Figure 3b). The resolutions we use in MAR have minimal effect on the timing and magnitude of the BIE. The initial modeled BIE is slightly lower than the observed BIE, indicating that the 0.6 threshold for band 2 in MOD09GA could be too high or that the onset of bare ice is delayed in MAR relative to the observed onset. The lower initial modeled BIE may also suggest that MAR 25 misses uncovering of ice by snowdrift which would blow snow into crevasses and holes. In early July, the modeled BIE quickly surpasses the minimum observed estimate. The late-July dip in modeled BIE is in part due to snowfall events over the BIE as simulated by MAR.  (Wilks, 2011).

Table 1: Contingency table of forecasts (MAR) and observations (MODIS) in terms of
determining either bare ice or snow. Bare ice is both forecast and observed is a hit (a), bare ice is forecast but snow is observed is a false alarm (b), snow is forecast but bare ice is 5

Bare ice albedo
The average observed bare ice albedo over 2000-2021 exhibits high spatial variability   shown to play a minor role in lowering albedo due to its sparse spatial distribution (Ryan et al., 2018). A part of the seasonal decrease in bare ice albedo also arises from accumulated surface meltwater on the bare ice, which may be misrepresented in MAR.

Meltwater production
20 Figure 6a and 6b show annual and seasonal averages, respectively, of the fractional contributions of each of the constituents in the meltwater production parameterization, where we use albedo values from MODIS in the shortwave radiation term. The shortwave radiation term is consistently the largest term contributing to meltwater production owing to the low albedo of bare ice and the long days during boreal summer. Shortwave radiation contributes on average 5 -5.5 25 times more to meltwater production than sensible heat flux, the second largest contributor.
Longwave radiation contributes to a net loss of heat in general during the study period, releasing more energy from the surface of the ice than it absorbs. This leads to a net negative contribution to meltwater production. Sensible heat flux contributes a small but rather constant fraction to meltwater production. The contribution of latent heat flux to meltwater production is very small and often negative (-0.11 -0.02), meaning that the latent heat flux on average results more in evaporation and sublimation than condensation and deposition. We quantify the effect of the bare ice albedo bias by determining the meltwater production in two scenarios: once with the observed albedo and once with the modeled albedo as input to the meltwater production parameterization. We use the ratio of the daily averages of these two meltwater production estimates to isolate the effect of bare ice albedo on meltwater production from the bare ice zone below 70° N on a daily basis (Figure 7). A positive ratio indicates that using 5 the modeled albedo in the parameterization results in an underestimation of the average meltwater production on that day, compared to using the observed albedo. We observe a strong increasing seasonal trend in the average seasonal ratio of 0.0245 per week from June through August with seasonal average ratios for June, July and August of 1.32, 1.42, and 1.54, respectively, indicating a seasonally increasing underestimation of parameterized meltwater production using the modeled albedo. The increase in the seasonal average ratio of 15 meltwater production indicates that the spatial distribution and intensity of LAC concentrations on bare ice could be increasing during the season, reaching peak values in late August. Increasing amounts of accumulated surface meltwater could also be a cause of the increasing trend in the ratio. As shown earlier, the daily bare ice extent significantly decreases after mid-August. This means there is a smaller area over which meltwater production is calculated, increasing its variability. The seasonal average reaches a low of 1.19 on June 7, and a peak of 1.75 on August 24. The meltwater production ratio exhibits significant daily variability throughout the study period. Daily ratios in June vary from 0.34 to 3.04, though these extreme ratios are sporadic. Since the overlapping bare ice extent between observation and model is small in early June, this period shows extreme variability in observed bare ice albedo and, thus, in meltwater production estimates. 5 Hence, strong conclusions cannot be drawn for early June. The minimum and maximum daily ratios in July are 0.81 and 2.11, respectively. August exhibits numerous extremely high daily ratios, especially in late August, with minimum and maximum ratios of 0.98 and 3.09, respectively.
The average annual ratio exhibits an increasing trend of 0.043 per decade from 2000 through 2021, indicating that parameterized meltwater production from the bare ice zone below 70° N is being 10 increasingly underestimated when using modeled albedo versus observed albedo. One explanation for this could be an annually increasing amount of LACs that are deposited onto and exposed in the bare ice zone during the summer. Increasing temperatures and accumulated surface meltwater could also create more favorable conditions for algal bloom growth. The minimum and maximum annual average ratios are 1.26 and 1.58, in 2004 and 2019, respectively. Averaged over the bare 15 ice zone below 70° N and the entire study period, the MAR-derived meltwater production using the modeled albedo could be underestimated by 42.8%, owing to an average overestimation of modeled bare ice albedo of 22.8%. The meltwater productions using observed and modeled albedo have a correlation coefficient of R 2 = 0.60.
In addition to the examination of time series, we average the meltwater production ratios 20 over the study period and map them onto the bare ice extent (Figure 8). We find high ratios over the dark ice zone in southwest Greenland, which is as expected from the high LAC concentrations in this area . Moreover, we find values between 0.4 and 1 higher up in the ablation zone in southwest Greenland. In this region, meltwater production estimates are close to 0 using both observed and modeled albedo. Melting also occurs significantly less frequently with 25 increasing elevation in this region. Hence, a small difference in albedo can result in a large percentage change of simulated meltwater production. We find extremely high ratios along the eastern margin where the observed albedo is significantly lower than the ~0.55 simulated by MAR ( Figure 4). The large albedo differences and low number of melting days in this area makes meltwater production estimates more variable.

5
An analysis by (Ryan et al., 2018) on sources of spatial albedo variability along a transect in the dark ice zone, perpendicular to the ice margin, found that 73% of spatial albedo variability can be attributed to LACs on the surface of bare ice; a mixture of algae, dust and black carbon.
Only 15% of the spatial variability of albedo is explained by accumulated surface meltwater in rivers, streams, ponds and lakes. Crevasses are responsible for 12% of the observed albedo 10 variability. Despite the very low albedo of cryoconite, due to its low abundance it accounts for only 0.6% of the albedo variability. Moreover, accumulated surface meltwater may act as a distributor of LACs; a small change in accumulated surface meltwater may thus result in larger albedo changes than merely the added surface water. Granted, the analysis by (Ryan et al., 2018) only holds for one transect covering 12.5 km 2 on August 6, 2014, and may not necessarily be representative of the entire bare ice area. A modeling study from (Goelles and Bøggild, 2017) suggests that melt-out of englacial black carbon and dust are responsible for most of the LACdriven meltwater production. Atmospheric deposition of black carbon and dust has a significantly lower effect on meltwater production in their study. These results hold for the location of an AWS (KAN_M) at 1270 m a.s.l, in the dark ice zone, for 2010-2015. It should be noted that biological 5 activity is not included in their model, so they could be interpreting some albedo changes due to biological activity as effects from black carbon and dust. It thus remains unclear how the effect of algae relates to the effect of dust and black carbon on bare ice albedo and meltwater production.
Though, qualitatively, we assume the conclusions drawn by (Goelles and Bøggild, 2017;Ryan et al., 2018) hold for the albedo differences between model and observations we find in our analysis. 10 We recognize that the choice of the threshold in band 2 of MOD09GA to determine the bare ice extent adds additional uncertainty to our results. (Shimada et al., 2016)  for this threshold are available in the current literature. We therefore believe that the 0.6 threshold is currently the best estimate.
We want to emphasize that our results on albedo and meltwater production only hold for the bare ice extent below 70° N that is simultaneously observed by MODIS and modeled by MAR.
Since the MOD10A1 product may be less reliable above 70° N, we exclude this region from this part of our analysis. This means that the albedo and meltwater production results shown in this study cannot be extrapolated to the northern half of the GrIS, though we expect the physical processes to be fairly similar over the entire GrIS. An improved understanding of errors in satellitederived albedo measurements, including additional high-quality in situ measurements, would be 5 useful for properly analyzing the effects of albedo on meltwater production above 70° N.
The ratios we mention in Section 3.3 pertain to the parameterized meltwater production using observed albedo and modeled albedo. No direct conclusions can thus be drawn on the performance of MAR in simulating meltwater production over the bare ice zone. However, in the Methods section we show that the parameterized meltwater production using the modeled albedo 10 and the original meltwater production in MAR have a very high correlation (R 2 = 0.92). We therefore believe that our conclusions are likely transferable and applicable to the performance of MAR in simulating meltwater production.
It is also important to note that the SMB simulated by MAR compares very well with SMB observations from the Programme for Monitoring the Greenland Ice Sheet (PROMICE) on average 15 over the entire GrIS . This suggests that the effects on meltwater production of a too high bare ice albedo through absorption of shortwave radiation might be compensated by other energy fluxes (LWnet, SHF, LHF) in the energy balance equation of MAR over the bare ice area. This is discussed in , who highlighted an overestimation of albedo and downward shortwave radiation but an underestimation of downward longwave radiation. 20 However, MAR underestimates melt at AWS locations in the ablation zone where melt is larger than 2 mWE/yr , suggesting that at these locations a lower bare ice albedo would improve comparison of modeled SMB with observations from PROMICE. At these locations, the SMB modeled by the Regional Atmospheric Climate Model (RACMO2.3p2) compares better with PROMICE than the SMB from MAR does, most likely since RACMO 25 integrates MODIS albedo into their model .

Conclusions
Using remote sensing observations and an energy balance model to parameterize meltwater production, we analyze the performance of the regional climate model MAR in simulating the spatio-temporal variability of the bare ice extent and albedo as observed by MODIS. We have shown that MAR performs reasonably well in simulating the bare ice extent on an annual basis.
Despite the similarities in maximum annual bare ice extent, MAR overestimates the daily bare ice extent during peak bare ice season from mid-July through mid-August. We also conclude that MAR overestimates bare ice albedo below 70° N on average by 22.8% during the study period. 5 This complements and builds further on a study by Alexander et al. 2014, who analyzed surface albedo over the entire GrIS. We advocate that this significant difference in bare ice albedo arises in substantial part from the lack of LAC representation in MAR's bare ice albedo scheme. A misrepresentation of accumulated surface meltwater on bare ice in MAR may also in part cause the difference between observed and modeled bare ice albedo. Using the meltwater production 10 parameterization, we isolate the effect of the bias in observed and modeled bare ice albedo on the meltwater production from the bare ice zone below 70° N. We find that, using the modeled albedo in the parameterization, meltwater production is underestimated on average by 42.8% during the study period. The underestimation of meltwater production increases on average with 2.45% per week from June through August and with 4.3% per decade from 2000 through 2021. The largest 15 discrepancies in meltwater production are located over the dark ice zone, where the highest LAC concentrations are found, and along the eastern margins of the ice sheet, where simulating bare ice extent is more difficult owing to the steep topography of the fjords and cliffs. Since meltwater production estimates from MAR and estimates from the parameterization with the modeled albedo are closely linked (R 2 = 0.92), we believe that the results pertaining to meltwater production are 20 likely transferable and applicable to MAR's performance in simulating meltwater production.
The results of this study show that research efforts should be directed towards uncovering the spatial and temporal variability of the distribution and trends of LAC concentrations on bare ice. Regional climate models, such as MAR, should work towards adopting a bare ice albedo scheme that allows for inputting spatially and temporally variable LAC concentrations on bare ice. 25 Radiative transfer models such as the SNow ICe and Aerosol Radiative model SNICAR are being improved to allow for inputting black carbon, brown carbon, dust, ash, and algae with a range of properties in a variable concentration (Whicker et al. 2021, in press). However, no GrIS-wide, or dark ice zone-wide, quantification of distributions and trends are available yet. Hence, this aspect of LACs on bare ice still has to be parameterized.
As global, and more so Arctic, atmospheric temperature continues to rise, more bare ice will be exposed by melting the snow that usually blankets the bare ice, increasing meltwater production from the ablation zone of the GrIS. Correctly modeling and predicting bare ice albedo, and in particular LAC concentrations on bare ice, is thus becoming increasingly imperative for proper projections of meltwater production from the GrIS by regional climate models and general 5 circulation models.

Code and data availability
The code underlying the processing and analysis described in this paper are available from Raf Antwerpen upon request. The MARv3.12 code and outputs are available at Author contributions 15 RA, MT, PA and WB designed the study. XF ran MARv3.12 simulations and provided the outputs. RA processed the data, performed the analysis, generated the figures and wrote the manuscript. All authors discussed the results and contributed to the final manuscript.

Competing interests
The authors declare they have no conflict of interest.