The Effect of Physically Based Ice Radiative Processes on Greenland Ice Sheet Albedo and Surface Mass Balance in E3SM

A significant portion of surface melt on the Greenland Ice Sheet (GrIS) is due to dark ice regions in the ablation zone, where solar absorption is influenced by the physical properties of the ice, light absorbing constituents (LACs), and the overlying crustal surface or melt ponds. Earth system models (ESMs) typically prescribe the albedo of ice surfaces as a constant value in the visible and near‐infrared spectral regions. This work advances ESM ice radiative transfer modeling by (a) incorporating a physically based radiative transfer model (SNow, ICe and Aerosol Radiation model Adding‐Doubling Version 4; SNICAR‐ADv4) into the Energy Exascale Earth System Model (E3SM), (b) determining spatially and temporally varying bare ice physical properties over the GrIS ablation zone from satellite observations to inform SNICAR‐ADv4, and (c) assessing the impacts on simulated GrIS albedo and surface mass balance associated with modeling of more realistic bare ice albedo. GrIS‐wide bare ice albedo in E3SMv2 is overestimated by ∼4% in the visible and ∼7% in the near‐infrared wavelengths compared to the Moderate Resolution Imaging Spectroradiometer. Our bare ice physical property retrieval method found that LACs, ice crustal surfaces, and melt ponds reduce visible albedo by 30% in the bare ice region of the GrIS ablation zone. The realistic bare ice albedo reduces surface mass balance by ∼145 Gt, or 0.4 mm of sea‐level equivalent between 2000 and 2021 compared to the default E3SM. This work highlights the importance of simulating bare ice albedo accurately and realistically to improve our ability to quantify changes in the GrIS surface mass and radiative energy budgets.


Introduction
Global mean sea level is expected to continue to rise in the coming century, and the Greenland Ice Sheet (GrIS) is currently the largest cryospheric contributor to increasing sea level (Casey et al., 2017;Fox-Kemper et al., 2021;Goelzer et al., 2020;King et al., 2020;McCutcheon et al., 2021;Steger et al., 2017).The total mass loss from the GrIS can be separated into two major categories: solid ice losses from discharge along the grounding line and surface runoff (Mouginot et al., 2019;van den Broeke et al., 2016).Since 2000, GrIS mass loss from surface runoff has exceeded that from solid ice discharge (King et al., 2020;Mouginot et al., 2019;van den Broeke et al., 2016).The majority of surface melt comes from the ablation zone, which constitutes a relatively small portion of the ice sheet (Box et al., 2012;Steger et al., 2017).Surface mass loss is projected to increase through the 21st century.However, the total projected mass loss from the GrIS is uncertain, partially due to our limited ability to accurately model the processes that contribute to the surface mass balance, or the annual difference between accumulation and ablation (Cogley et al., 2010;Goelzer et al., 2020;Otosaka et al., 2023;van den Broeke et al., 2008).
Ice melt on the surface of the GrIS is in part modulated by the surface albedo.High albedo surfaces, such as fresh snow, efficiently reflect insolation, while dark surfaces, such as dirty bare ice, absorb a large portion of the incoming shortwave radiation.In the ablation zone of the GrIS, the albedo and insolation strongly control the surface mass loss during the summer (Tedesco et al., 2011;van den Broeke et al., 2008).The ablation zone has a lower albedo than the interior of the ice sheet due to surface melt, retreat of the snowline elevation exposing darker bare ice, and contamination by light absorbing constituents (LACs), such as black carbon, dust, and glacier algae (Cook et al., 2020;Ryan et al., 2018Ryan et al., , 2019;;Tedstone et al., 2017Tedstone et al., , 2020;;Williamson et al., 2018).A significant portion of surface melt has been attributed to the growth of glacier algae on the ice sheet (Cook et al., 2020;Stibal et al., 2017;Wang et al., 2018;Williamson et al., 2018Williamson et al., , 2020;;Yallop et al., 2012).These algae present a positive feedback, as they darken the surface of the ice sheet with their photoprotective pigments, leading to increased surface melt, releasing nutrients frozen in the ice, and allowing more glacier algae to grow (Di Mauro et al., 2020;McCutcheon et al., 2021;Williamson et al., 2019).
Significant melt has also been attributed to the exposure of bare ice (Antwerpen et al., 2022;Ryan et al., 2019;Tedstone et al., 2017).The snowline elevation feedback, or the seasonal exposure of bare ice as the snowline retreats up the ice sheet during the melt season, has been found to amplify ice sheet melt up to five times more than hydrological and biological processes between 2001 and 2017 (Ryan et al., 2019).This is because bare ice surfaces are generally much darker than snow surfaces, with bare ice albedo generally ranging from 0.1 to 0.6 (Barry, 1996;Whicker et al., 2022;Wiscombe & Warren, 1980).This increases the absorption of solar radiation, which leads to increased surface melt.Ice surfaces are also much less porous than snow surfaces, therefore, melted runoff is much less likely to be captured and refrozen once the porous surface snow and crust have melted off.
Over the past few decades, the GrIS has darkened and is expected to continue to darken with increasing polar temperatures (Shimada et al., 2016;Tedesco et al., 2016).The spatial and temporal extent of the ablation zone has already increased due to higher snowline elevations and more expansive glacier algae colonies (Colosio et al., 2021;Shimada et al., 2016;Tedesco et al., 2016).The darkening of the ice sheet and lengthening of the melt season both present positive feedbacks, propagating more surface melt.Therefore, it is increasingly important to accurately model the influence of bare ice and LACs on GrIS albedo.
Significant progress has been made on our ability to model bare ice surfaces and estimate the total impact of bare ice exposure on GrIS melt.Previous work has explicitly represented bare ice using the physical microstructure of the ice, including the size distribution of air bubbles within the ice, the density of the ice, and a refractive boundary that represents the transition in the index of refraction between air media and ice media (Briegleb & Light, 2007;Dadic et al., 2013;Gardner & Sharp, 2010;Mullen & Warren, 1988;Robledano et al., 2023;Whicker et al., 2022).Subsequent work further developed an offline snow radiative transfer model to accurately simulate a spectrally resolved cryospheric column of snow and ice with a refractive boundary and incorporates LACs, including black carbon, and algae (Whicker et al., 2022).This modeling work provides a useful theoretical range of bare ice albedo and compares well to in-situ ice measurements.However, the new ice radiative transfer techniques have not been applied to estimate ice sheet-wide impacts on albedo associated with bare ice exposure.
Satellite data have been used to estimate ice-sheet-wide albedo and melt trends.Antwerpen et al. (2022) compare the regional climate model Modèle Atmosphérique Régional (MAR) with satellite imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS) for the GrIS south of 70°N and found that MAR potentially overestimates bare ice albedo by 22.8%, which leads to a 42.8% underestimation of surface melt in the bare ice zones between 2000 and 2021.Other studies have used MODIS albedo observations to prescribe bare ice albedo on the GrIS in regional models, such as the Regional Atmospheric Climate Model (RACMO2.3p2), the HIRHAM model (Antwerpen et al., 2022;Box et al., 2012;Langen et al., 2017;Noël et al., 2018Noël et al., , 2019)).Satellite data has also been used to infer information about snow surfaces, most typically the specific surface area or effective snow grain radius.The physical properties of snow and ice are typically inferred using reflectance measurements in the near-infrared (NIR) region of the spectrum, which is sensitive to changes in grain size and relatively insensitive to the influence of LACs (Dozier et al., 1981(Dozier et al., , 2007;;Flanner et al., 2021;Painter et al., 2012).New work is utilizing satellite imagery and imaging spectroscopy to infer snow and ice properties such as the specific surface area, LACs content, and even snow and glacier algae concentrations (Bohn et al., 2021(Bohn et al., , 2022;;Cook et al., 2020;Kokhanovsky et al., 2019;Painter et al., 2001Painter et al., , 2012;;Vandecrux et al., 2022;Wang et al., 2018Wang et al., , 2020)).Other studies have used satellite measurements to quantify spatial and temporal changes in the ablation zone of the ice sheet (Alexander et al., 2014;Antwerpen et al., 2022;Navari et al., 2021;Shimada et al., 2016;Tedesco et al., 2013Tedesco et al., , 2016;;Tedstone et al., 2017).For example, MODIS products have been used to categorize the spatial extent of bare ice on the GrIS and to compare to regional climate models' albedo (Alexander et al., 2014;Antwerpen et al., 2022;Ryan et al., 2019;Shimada et al., 2016;Tedstone et al., 2017).Significant scientific progress has been made in both the modeling and measuring of ice processes, however, global Earth System Models (ESMs) such as the Community Earth System Model (CESM) and the Energy Exascale Earth System Model (E3SM) both prescribe constant bare ice albedo values in the visible and near-infrared (NIR) regions of the spectrum, thereby neglecting spatial and temporal albedo variability shown in observations.This work improves the current representation of bare ice albedo on the GrIS within E3SM by implementing the bare ice radiative transfer modeling technique developed by Whicker et al. (2022) within the E3SM land model (ELM).This work utilizes MODIS bare ice albedo measurements and applies snow and ice remote sensing techniques to MODIS data to infer properties about the specific surface area and LAC content of the regions categorized as bare ice surfaces.The spatially and temporarily varying MODIS inferred bare ice properties are utilized within SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) in ELM, unlike the implementation of SNICAR in the land snow model, which simulates evolution of snow grain size.ELM is forced with a data atmosphere to calculate the GrIS ice albedo and corresponding surface mass balance from 2000 to 2021.Our results indicate the default ELM bare ice albedo values overestimate the exposed bare ice albedo and underestimate total surface melt from the GrIS compared to the spatio-temporally varying bare ice albedo simulations, particularly in the southwest ablation zone, due to neglecting realistic varying bare ice albedo.These are the first ESM simulations that contain improved ice radiative transfer physics that allow for spatial and temporal variation in ice albedo.It is also the first empirically constrained and physically based estimate of the bare ice and impurity effects on the GrIS ablation zone radiative and melt budget.This work also introduces a novel ice physical property retrieval algorithm used to inform the ESM.
The Models, Data, and Methods section provides background on the models and data utilized in this work.It also outlines how the SNICAR-ADv4 radiative transfer model is implemented within ELM and how various MODIS products are used to inform SNICAR-ADv4 bare-ice albedo calculations.The Results section compares the MODIS-informed SNICAR-ADv4 enabled simulation to the default E3SM simulation and discusses the impact of spatially and temporally varying bare ice albedo on the radiative budget and the surface mass balance.Overall, this work highlights the growing importance of accurately simulating the albedo of exposed bare ice on the GrIS and provides novel techniques for modeling and quantifying the impact bare ice albedo has on regional energy and mass balances.

Models: The SNICAR-ADv4 Model
The SNow, ICe, and Aerosol Radiation model with the Adding-Doubling technique Version 4 (SNICAR-ADv4) is a multilayer two-stream radiative transfer model for cryospheric surfaces that uses 480 wavelength bands between 0.2 and 5 μm.Dang et al. (2019) revised the previous snow radiative transfer scheme within ELM that utilized the snow optical properties derived from Wiscombe and Warren (1980), and Warren and Wiscombe (1980) and the two-stream solution derived by Toon et al. (1989) to instead utilize the delta-Eddington adding-doubling two-stream approximation presented in Briegleb and Light (2007).Dang et al. (2019) found that the delta-Eddington adding-doubling two-stream approximation results in the most accurate albedo compared to other two-stream methods and allows for the incorporation of internal Fresnel layers (Briegleb & Light, 2007).The ability to incorporate a refractive boundary using the Fresnel equations enables simulating columns with snow, liquid water, and solid ice (Briegleb & Light, 2007;Whicker et al., 2022).Whicker et al. (2022) developed SNICAR-ADv4 by extending Dang et al. (2019)'s snow albedo model to simulate the albedo of snow and ice by including Fresnel reflection between air and ice media (Briegleb & Light, 2007) and scattering air bubbles within an absorbing ice medium (Gardner & Sharp, 2010;Mullen & Warren, 1988) (Figure 1).The use of SNICAR-ADv4 allows for the calculation of bare ice albedo using the physical properties of the ice, such as the specific surface area (SSA), the ice density (or air volume fraction), and effective air bubble radii within the ice, and the atmospheric conditions (Whicker et al., 2022).The SSA is defined as the ratio of the total surface area of ice-air interfaces to its mass and is directly related to snow and ice albedo (Whicker et al., 2022).In SNICAR-ADv4, snow layers are represented as ice crystals within air, and ice layers are represented as solid ice layers with air bubble inclusions.SNICAR-ADv4 requires layer-specific information about the snow grain size, shape, and mass path, the air bubble radii, and ice density.SNICAR-ADv4 also includes LACs such as dust, black carbon, snow algae, and glacier algae (Chevrollier et al., 2022;Cook et al., 2020;Flanner et al., 2007Flanner et al., , 2021;;Skiles et al., 2018;Whicker et al., 2022).This allows the model to quantify the impact of LACs on spectral albedo.

Models: The Energy Exascale Earth System Model (E3SM)
The Energy Exascale Earth System Model version 2 (E3SMv2) is the Department of Energy's ESM (Golaz et al., 2022).The land component of the E3SM model, ELM, was based on the Community Land Model version 4.5.The default ELM snow and ice layering scheme over land ice is described by Oleson et al. (2013).Sea ice uses a different, but related radiative transfer scheme that is outside the scope of the current study (described in Briegleb & Light, 2007).The snow radiative transfer scheme (SNICAR-AD) uses the physical properties of the snowpack and albedo of the top layer of the underlying ground to determine the column albedo (Dang et al., 2019;Flanner & Zender, 2006).ELM contains an advanced model for determining snow properties from incoming precipitation and snow grain growth theory (Flanner & Zender, 2006;Hao et al., 2022).In the case of land ice, the underlying albedo is a constant 0.6 in the visible and 0.4 in the NIR (Oleson et al., 2013).The constant albedo values are propagated into the calculations of fluxes and ice melt.

Methods: Implementing SNICAR-ADv4 Into ELM
This work incorporates SNICAR-ADv4 into the land component (E3SM Land Model, ELM) of the Energy Exascale Earth System Model version 2.1 (E3SM) (Figure 1).Utilizing SNICAR-ADv4 within ELM allows for the radiative processes to be simulated over ice surfaces as well as snow surfaces.To implement SNICAR-ADv4 into ELM, the offline 480 band spectral resolution of the bare ice optical properties (air bubbles within ice, and LACs) are weighted to the SNICAR in ELM five-band spectral scheme using incoming spectral flux simulated with a typical sub-Arctic summer clear sky atmospheric profile.We also introduce the Fresnel equations into SNICAR-AD in ELM to simulate the changing index of refraction between snow surfaces and ice surfaces.We extend the SNICAR-AD radiative transfer model within ELM to represent glacier ice by utilizing the vertical ice layers within ELM land ice grid cells.These new radiative transfer physics allow light to penetrate the ice surface and calculate ice albedo based on the ice physical properties and environmental conditions.However, the solid ice layers in land ice grid cells do not yet include properties that are relevant for our radiative transfer calculations (e.g., density or LAC content).The physical properties of ice are either unknown (e.g., air bubble radius and SSA) or are treated as constant (ice density) within ELM, thus we utilize satellite observations to estimate the physical properties of bare ice on the GrIS.

Data: MODIS Reflectance and Albedo Products
MODIS reflectance (MOD09GA) and albedo (MCD43A3) products are used from the summer melt seasons (June, July, and August (JJA)) from 2000 to 2021 to retrieve the optically relevant physical properties of bare ice surfaces.The MCD43A3 products are composed of atmospherically corrected MODIS reflection measurements collected over 16-day moving average periods to derive daily integrated albedo at local solar noon.While the Black Sky Albedo (BSA) and White Sky Albedo (WSA) are both available, the calculation of blue sky albedo requires a model or observations of the atmospheric state (Stroeve et al., 2005).Previous work has shown there is little difference between BSA and WSA for a typical summer day and the use of BSA is acceptable for analyses over the GrIS during JJA (Alexander et al., 2014;Stroeve et al., 2005).Thus, this work utilizes the MCD43A3 BSA band 2, 8, and 9 albedo products.

Methods: Upscaling and Filtering the MODIS Products
For the purpose of this study, the MODIS products were upscaled from their 500m resolution to the ELM 0.5 × 0.5-degree resolution grid using a method that bins each satellite pixel contained within a given model grid and averages the pixel values, weighting each image equally.The MODIS-supplied cloud mask from MOD10A1 was used to omit pixels with cloud cover and a bare ice mask was used to omit pixels that were determined to be snow-covered (Antwerpen et al., 2022).Following the methodology presented by Shimada et al. (2016) and utilized by both Antwerpen et al. (2022) and Tedstone et al. (2017), the MOD09GA band 2 (841-876 nm) reflectance product is used to determine where bare ice is exposed on the ice sheet.These studies develop a satellite-derived bare-ice extent by setting a threshold of 0.6 to the MOD09GA band 2 reflectance product.Any pixels with a reflectance less than or equal to 0.6 are defined as bare ice.It is important to note the uncertainty associated with this thresholding technique is unknown, however, a comparison with Landsat imagery shows good agreement and a relative error of 0.16%, with the surface classification (Antwerpen et al., 2022).The bare ice extent is also filtered to exclude pixels with surface elevations higher than the mean equilibrium line altitude of 1679m a.s.l. for 2000-2021, as determined by Antwerpen et al. (2022).While the total number of valid bare ice MODIS pixels varies based on snow cover and data availability, the output 0.5 × 0.5-degree daily grid cell albedo averages contain all the available bare ice MODIS measurements within a given grid cell.To inform the SNICAR-ADv4 ice properties within ELM we utilize MODIS observations and relationships between albedo, SSA, and the physical properties of ice (e.g., Figure 2, described next).

Methods: Inferring Ice Physical Properties From MODIS Bare Ice Albedo
The 480-band offline SNICAR-ADv4 is used to generate a lookup table of bare ice physical properties and the corresponding SSA and albedo.The lookup table was generated by running the offline SNICAR-ADv4 over varying bare ice conditions.The bare ice physical properties utilized within the lookup table are constrained using in-situ measurements of ice properties (Chevrollier et al., 2022;Cooper et al., 2018;Cuffey & Paterson, 2010;Dadic et al., 2013).The bare ice physical properties in the lookup table simulations are restricted to ice densities ranging from 650-916 kg/m3 and air bubble radii ranging from 100-1500 μm.The offline simulations for the lookup table were run over a range of SZAs (35°-75°) that correspond to ELM GrIS grid cells' local noon SZAs to be consistent with the MCD43A3 albedo products.The range in the SNICAR-ADv4 offline lookup table band 2 albedo as a function of SSA is displayed in Figure 2.
After the daily MODIS albedo products are filtered and upscaled, an inverse modeling technique is employed using the MCD43A3 band 2 albedo and the offline SNICAR-ADv4 lookup table to find a set of ice physical properties that most closely match MODIS albedo.The inverse model takes the MCD43A3 band 2 BSA product (which is in the NIR part of the spectrum), and the grid cell local noon SZA as the initial inputs.The inverse model compares each grid cell averaged MCD43A3 band 2 BSA to the lookup table band 2 albedo from all possible bare ice properties with the same local noon SZA as the upscaled MODIS grid cell.Once the closest albedo value is found in the lookup table for a given MCD43A3 band 2 BSA and SZA, the corresponding bare ice physical properties (SSA, ice density, and air bubble radius (R eff )) are retrieved from the lookup table for that albedo.This results in a SSA, ice density, and effective air bubble radius, for the MCD43A3 band 2 albedo values in each grid cell from the offline lookup table, as shown in Figure 3.These physical properties are then fed into SNICAR-ADv4 in ELM.
It is important to note, the physical properties of bare ice can be determined from the SSA, however, the relationship between SSA and ice density and air bubble radius is not unique (Whicker et al., 2022: equation 5).There is essentially a one-to-one relationship between SSA and albedo for a given SZA, meaning the combination of ice density and bubble size selected to achieve a desired SSA is of little practical importance for the radiative transfer modeling (Figure 2).Thus, whichever set of ice physical properties result in the closest band 2 agreement are selected.It is also worth noting, that the NIR part of the spectrum is used to determine the physical properties of bare ice because it is sensitive to changes in snow and ice morphology and is minimally influenced by LACs that are present on the ice surface (Dozier et al., 1981(Dozier et al., , 2007;;Painter et al., 2012;Schneider et al., 2019).This allows for the physical properties of the ice to be inferred without the influence of darkening LACs.
The resulting map of retrieved bare ice properties accurately reproduces (within 2.9%) the MODIS NIR broadband (0.7-5.0 μm) albedo (Figure 4).However, these bare ice properties significantly overestimate the bare ice albedo in the visible part of the spectrum in the offline lookup table compared to the MCD43A3 visible albedo (Figure 5).This is expected because the visible part of the spectrum is strongly influenced by LACs, whereas the NIR spectrum is not.The inverse modeling technique initially utilizes the NIR part of the spectrum to inform the SSA (and ice density and air bubble effective radius).However, this does not provide any insight into the visible part of the spectrum, where the influence of LACs, melt ponds, and crustal surfaces dominate the albedo.
To account for the influence of LACs on bare ice and achieve a visible albedo similar to the MCD43A3 product, we apply the upscaling and inverse modeling technique to the visible MCD43A3 albedo and assign an equivalent black carbon concentration from the offline lookup table that resolves the discrepancy between the SNICAR-ADv4 inverse modeled albedo and the MCD43A3 visible albedo.The equivalent black carbon concentration is determined using the ice properties previously obtained from the NIR region of the spectrum and then finding the concentration of black carbon in ice with the corresponding SSA that minimizes the total difference between visible MCD43A3 albedo and offline SNICAR-ADv4 lookup table visible albedo.
Obtaining the SSA of the ice before deriving the LAC concentration is important because the efficacy of darkening agents depends on the SSA of the ice.The difference in visible albedo between clean ice with the NIRdetermined ice physical properties and the ice that contains the equivalent black carbon with the same ice physical properties is the "Impurity Effect" which is described in more detail in the results.In reality, a variety of different LACs and surface conditions contribute to the darkening of bare ice (Bøggild et al., 2010; Chevrollier et al., 2022; Dadic et al., 2013;Tedstone et al., 2020;Whicker et al., 2022), but our primary goal here is merely to derive a set of ice properties that, when prescribed in ELM, result in an improved representation of the spatial and temporal variability of bare ice albedo.
From the MODIS data and the offline SNICAR-ADv4 lookup table, we produce maps at the ELM r05 grid scale (0.5°× 0.5°-degree) of (a) specific surface area, (b) bare ice density, (c) air bubble radii, and (d) equivalent black  carbon concentration that vary spatially over the GrIS and temporally (monthly in JJA from 2000 to 2021) (Figures 3 and 5).A set of fill parameters are specified over any regions for which there is no MODIS bare ice data, that is, where there was not at least one valid bare ice pixel for each summer month.To minimize the impact of the fill albedo values on the analysis below, the fill bare ice properties are assigned based on ice properties that yield the closest albedo to the default E3SM constant visible albedo of 0.6 and NIR albedo of 0.4 in the offline lookup table at a SZA of 73°, which corresponds to the ELM GrIS area and cosine weighted annual average SZA.With the specification of ice properties rather than ice albedo, this allows the bare ice albedo to vary more realistically with SZA and atmospheric conditions instead of being held constant at 0.6 in the visible and 0.4 in the NIR.It is worth noting that the majority, if not all, of these fill pixels, are in snow-covered glacier areas, for example, above the equilibrium line, where the ice albedo has minimal or no impact on the surface albedo.These data sets are then merged into an input file ingested by ELM.

Methods: Implementing the MODIS Informed Bare Ice Properties in SNICAR-ADv4 Enabled ELM Simulations
The monthly MODIS-informed bare ice properties are applied in ELM during May-September, ensuring that the bare ice properties are properly weighted during June-August in the ELM time interpolation scheme.This is because SNICAR-ADv4 within ELM uses a mid-month temporal weighting scheme.For example, on June 15 the model uses the June ice properties (recorded as June 15), but on June 1, it is using a weight that is ½ the May value and ½ the June value, as June 1 is halfway between May 15 and June 15.The equivalent BC concentration only influences the albedo in the visible region of the spectrum to avoid any contamination in the NIR region, where ice properties were tuned to match MODIS observations without LAC influence.This constraint is quite minor because the influence of LACs on NIR albedo is generally very small, especially in the case of BC.Because SNICAR-ADv4 within ELM only uses 5 spectral bands, instead of the 480 bands the offline model uses, and the rest of the E3SM model only uses two bands (one visible and one NIR) a tuning parameter is introduced in the weighting of the 4 SNICAR-ADv4 NIR spectral bands into one broadband NIR albedo to ensure good agreement between the SNICAR-ADv4 ELM and MODIS bare ice albedo.This tuning parameter was determined using sensitivity test simulations over a random 2-year period within the study period .For any land icecovered ELM grid cell, two albedo values are calculated.First, the albedo of a bare ice column is calculated.Then, that albedo is used to replace the constant bare ice albedo value previously employed in ELM (0.6/0.4).Then, if snow is present over the bare ice, the albedo of the snow is calculated, using the new bare ice albedo as the lower boundary beneath the snow column.The total grid cell albedo is then calculated based on the fractional coverage of the different land types and snow cover.
The simulations are run with E3SMv2 with a 16-layer firn model in a land-only configuration with MPAS-Albany Land Ice model (MALI) on a 0.5 × 0.5-degree resolution and 10 elevation classes over the ice sheets (to account for variable elevation and therefore temperatures within a single model grid cell).For the first time, the European Center for Medium-Range Weather Forecast's fifth-generation atmospheric Reanalysis (ERA5) was used to derive a 6-hourly atmospheric forcing data set applied as the land surface boundary conditions to ELM (Hersbach et al., 2020).ERA5 was used because its precipitation rates over the GrIS are more highly correlated with measured net accumulation rates than other ELM atmospheric forcing options during the 1980s, a proxy period for the preindustrial Greenland surface climate (Schneider et al., 2023).Improved precipitation rates coupled with a more realistic firn density parameterization enable a more reliable calculation of the GrIS surface and nearsurface mass and energy budgets in simulations that determine ELM's Greenland's bare ice regions (Schneider et al., 2022;van Kampenhout et al., 2017).Historical simulations, starting in 1990, were initialized using conditions after a 300-year spin-up run, to reach a stable and realistic snowpack depth on the GrIS, for which the surface boundary conditions were provided by the 1980-1989 ERA5 forcing data set (repeating every 10 years).From 1990 to 1997, the simulation used the default ELM ice properties (ice albedo of 0.6 in VIS and 0.4 in NIR, e.g., see left side of Figure 1).In our experimental runs from 1998 to 2000, the 2000 JJA MODIS-informed bare ice properties were applied to a short spin-up run for the varying bare ice to avoid any strange model behavior between switching bare ice schemes.From 2000 to 2021 the spatially and temporally varying bare ice albedo was calculated in ELM by the MODIS-informed SNICAR-ADv4.For the control simulation, the 1990s simulations continued through 2021 with the default constant ice albedo from 1980 to 2021 and the same ERA5 data atmosphere forcing for the corresponding time period.

Journal of Geophysical Research: Atmospheres
10.1029/2023JD040241 Figure 6.Moderate Resolution Imaging Spectroradiometer (MODIS) albedo compared to the default ELM scheme and the SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) scheme.The albedo comparisons of (a-e) shortwave (0.3-5 μm) broadband albedo (BBA), (f-j) visible (0.3-0.7 μm), and k-o) NIR (0.7-5μm).From left to right in each row: the MODIS albedo (a,f,k), default ELM albedo (b,g,l), SNICAR-ADv4 enabled ELM albedo (c,h,m), the difference between the ELM default and MODIS albedo (d,i,n), and the difference between the SNICAR-ADv4 enabled ELM and MODIS albedo (e,j,o).In the difference plots, red values indicate that the ELM simulated albedo is larger and blue values indicate that the MODIS albedo is larger.The SNICAR-ADv4 enabled ELM simulations reproduce the MODIS albedo measurements accurately (within 0.007 of the MODIS albedo).
The results from 2000 to 2021 from the default ELM bare ice albedo simulation are compared to the results from 2000 to 2021 from the ELM SNICAR-ADv4 bare ice albedo simulation, and the simulated bare ice albedo from both simulations are compared with MODIS bare ice albedo.The ELM local noon bare ice albedo was calculated using the local noon downwelling and upwelling fluxes over bare ice.To isolate the influence of bare ice albedo, in the case of spatial analysis, the albedo of the snow above the ice is ignored, and in the case of spatially and temporally averaged numbers, the bare ice albedo is weighted by the bare ice area fraction.The GrIS surface mass balance (SMB) was calculated for both the experiment and control simulations to analyze the differences in the SMB between the two.

Results and Discussion
The SNICAR-ADv4 enabled ELM simulation utilizes bare ice physical properties that were inferred from MODIS products over the bare ice region of the ice sheet.These results present the first simulations using physically based radiative transfer methods over the GrIS.This section first describes the optically relevant bare ice physical properties inferred from MODIS data, then discusses the role LACs play in darkening the ablation zone.We present the albedo and surface mass balance results from our comparisons of the SNICAR-ADv4 enabled simulation and the default ELM simulation.

Bare Ice Properties
The MODIS albedo bare ice property retrieval reveals highly variable bare ice physical properties in the GrIS ablation zone (Figure 3).The SSA in the lookup table varies from 0.0024 to 26 m 2 kg 1 , while the ice density and effective radius (R eff ) vary within the ranges specified within the lookup table (650-916 kg/m 3 and 50-1500 μm respectively).Varying the ice physical properties that are applied in E3SM radiative transfer calculations results in much more realistic variability in simulated albedo (Figure 6).The multi-year average June, July, and August properties shown in Figure 3 could be applied as climatologies for ELM simulations of periods outside the MODIS era, similar to how climatologies of other land forcing data are applied.While the representativeness of these ice properties to time periods outside of the satellite era is unknown, the improved model physics enable a physically based treatment of ice radiative transfer processes and allow for the albedo of bare ice to change with changing ice physical properties and atmospheric conditions such as insolation and solar zenith angle (SZA).

Impurity Effect
The impurity effect is defined as the difference in the visible albedo determined by the NIR SSA with no impurities and the actual MODIS measured bare ice albedo.On the GrIS ablation zone, the bare ice area-weighted average reduction in visible albedo due to LACs is 0.3 (Figure 5).This estimation can serve as an initial upper bound for the influence of impurities and darkening surface structures on bare ice, as it simply measures the difference between the theoretical visible albedo determined from the ice physical properties and the actual measured albedo.Visible darkening may also result from surface roughness features that are unresolved in plane-parallel models like SNICAR (e.g., Larue et al., 2020) as well as melt ponds.The strength of the impurity effect varies along the ablation zone because it depends on the accumulation of LACs, the colonization of darkly pigmented glacier algae, and the formation of crustal surfaces and melt ponds.This effect is strongest along the terminal edges of the ice sheet.This is likely because the ice in these regions is exposed for longer durations of the summer melt season and is therefore able to accumulate more LACs both through deposition and the outcropping of dust and black carbon as ice melts, and has time to develop crustal surfaces and melt ponds (Gardner & Sharp, 2010;Warren, 1984; Warren &  Journal of Geophysical Research: Atmospheres 10.1029/2023JD040241 Wiscombe, 1980;Wientjes et al., 2012;Wientjes & Oerlemans, 2010).In the southwest ablation zone, it is likely that much of the darkening is caused by glacier algae (Bohn et al., 2022;Chevrollier et al., 2022;Cook et al., 2020;Ryan et al., 2018;Tedstone et al., 2017Tedstone et al., , 2020)).The influence of LACs, surface roughness, and melt ponds on albedo are all significant.Future work should attempt regional quantifications or constraints of the impact of each surface darkening agent on the entire GrIS ablation zone.

Summertime Albedo
In the summer melt season (June-August) the default prescribed ELM ice albedo is higher than the bare ice area weighted SNICAR-ADv4 enabled ELM bare ice albedo.The spatially and temporally average MODIS visible albedo is 0.451, 0.05 lower than that in the default ELM case (0.5) (Table 1).The SNICAR-ADv4 enabled simulation yields an exposed bare ice area-weighted average albedo of 0.450, with a difference of 0.009 from the MODIS value.These exposed bare ice-weighted albedo values are more typical of dark ice in the ablation zone.The uncertainty associated with the MCD43A3 product ranges from 0.02 to 0.05, with higher errors occurring over regions with complex terrain, and potential biases due to the data being limited to clear sky conditions (Campagnolo et al., 2016;Cescatti et al., 2012;Klok et al., 2003;Song et al., 2019;Stroeve et al., 1997;Wang et al., 2010Wang et al., , 2012;;Wen et al., 2022).This uncertainty is inherently propagated into the SNICAR-ADv4 enabled simulations, as the SNICAR-ADv4 enabled ELM bare ice albedo is tuned to achieve good agreement with the MCD43A3 albedo products used in the inverse modeling technique.
The SNICAR-ADv4 enabled simulations accurately reproduce MODIS albedo, as expected given that we derive bare ice properties from MODIS narrowband reflectance.Visible, NIR, and solar broadband albedo are all within 0.007 of the MODIS determined bare ice albedo (Figure 6).Some of the differences can be attributed to the coarse resolution spectral grid used in ELM and the weighting of the five ELM bands into two atmospheric bands, the visible (0.3-0.7 μm) and NIR (0.7-5.0 μm), which have been found to introduce errors of varying magnitudes (Wang et al., 2022).
Spatially, the default ELM bare ice values overestimate the albedo along the outer edge of the ice sheet, as this is where much of the dirty exposed bare ice is found but underestimates the albedo along the internal edge of the ablation zone (Figure 6).The default ELM's underestimation of bare ice albedo at higher elevations on the ice sheet can likely be attributed to the ice being exposed less frequently at higher elevations.When the higherelevation ice is exposed, it likely has a complex crustal surface and contains less impurities than ice at lower elevations that is exposed for more of the melt season (Tedstone et al., 2020).However, in ELM these regions are rarely, if ever exposed because the thick overlying snowpack translates into a high snow cover fraction (Figure 7).
There are no significant GrIS-wide temporal trends in exposed bare ice albedo (Figure 8), though bare ice albedo varies realistically with year in the experiment simulation, while not varying at all in the control simulation.The mean GrIS exposed bare ice albedo is consistently lower than the default ELM values, with the lone exception being June albedo in the visible spectrum (Figure 8).Bare ice albedo is overestimated in all other months by the default ELM scheme, especially in the NIR region of the spectrum, which is overestimated by at least 0.04 in 15 out of the 21 years simulated in this study (Figure 8).The average exposed bare ice BBA is highest in June, at the start of the melt season and continues to decline through July and August.However, the NIR albedo is similar in July and August, thus, the reduction in BBA between July and August is likely caused by increasing LACs and surface melt causing larger reductions in the visible part of the spectrum.

Surface Mass Balance
To assess the impact of the new SNICAR-ADv4 bare ice albedo scheme in ELM, we compare the SMB in the experiment and control simulations.For these comparisons, we utilize the Regional Atmosphere Climate Model (RACMO) GrIS mask because it is commonly used for GrIS SMB analyses and is similar to the ELM GrIS land ice region.It is worth noting the SMB calculations are sensitive to the ice mask used (Vernon et al., 2013).In the default ELM simulation, the 2000-2021 SMB is 368.7 Gt yr 1 with a (spatially averaged) standard deviation of 184 Gt (Figure 9), while in the SNICAR-ADv4 enabled simulations it is 362.1 Gt yr 1 (with a standard deviation of 189 Gt) meaning the MODIS informed varying bare ice albedo reduces the annual SMB by around 6.6 Gt yr 1 (∼2%).Over the course of the entire 21-year simulation, SMB is reduced by ∼145 Gt, or ∼0.4 mm of sea-level rise equivalent, in the SNICAR-ADv4 simulations.It is not accurate to assume all surface melt will immediately contribute to sea-level rise due to refreezing, percolation of meltwater into the ice sheet, and the formation of cryoconite holes.However, surface melt on bare ice surfaces is much more likely to contribute to sea-level rise, as ice is much less porous than snow and is able to capture much less liquid water than snow or firn surfaces (Ryan et al., 2019).
The default ELM bare ice albedo parameterization underestimates total melt from the GrIS compared to the MODIS informed SNICAR-ADv4 enabled simulations.These SMB estimates fall within the range of other studies that analyze the GrIS SMB over similar time periods; ranging from 456 to 529 Gt yr 1 for 1900-2015 reconstructions and 310-430 Gt yr 1 for 1990-2010 reconstruction, with a heavy dependence on the model, atmospheric forcing, and resolution used (Fettweis et al., 2017(Fettweis et al., , 2020;;Franco et al., 2012).However, the impact of the new albedo scheme varies spatially.Along the edge of the southwest ablation zone, the default ELM simulation underestimates mass loss (the blue region of the difference plot in Figure 9c) and along the inner edge of the southwest ablation zone, the default ELM simulation overestimates mass loss (the red region of the difference plot in Figure 9c).These results are consistent with the spatial trends in the bare ice albedo because the bare ice albedo values are buried beneath the snowpack and the higher albedo could be contributing to less total melt along the inner edge of the ablation zone in the SNICAR-ADv4 enabled simulation.
Just over 5 Gt yr 1 of the difference in SMB can be attributed to the summer melt season (JJA) (Figure 10).The total JJA SMB in the control simulation is 481.66Gt yr 1 , where the ELM bare ice albedo is constant.In the case of the experimental SNICAR-ADv4 simulation, the total JJA SMB is 486.72 Gt yr 1 .The average 1.53 Gt yr 1 difference between the annual average difference (6.59 Gt yr 1 ) and the JJA difference (5.06 Gt yr 1 ) can be attributed to the impacts of the new bare ice albedo scheme, as that is the only difference between the two simulations.The new albedo scheme reduces the SMB through the end of the melt season through September due to the implementation of the MODIS informed bare ice properties in May-September, where the MODIS inferred bare ice properties from June are utilized in May, and August values are utilized in September.This seasonal SMB cycle is similar to that found in previous work (Steger et al., 2017).The JJA spatial trends in the difference between the control and experiment simulations are similar to that in Figure 9, where we see an underestimation in mass loss along the outer edge of the southwest ablation zone, and an overestimation in mass loss along the inner edge of the southwest ablation zone in the control simulation (Figure 10).

Conclusions
An evaluation of the SNICAR-ADv4 enabled ELM simulation compared to the default simulation from 2000 to 2021 reveals that the default ELM method overestimates bare ice albedo by ∼5% in the shortwave broadband albedo (BBA) and underestimates total surface melt by ∼0.4 mm of sea level equivalent.These results highlight the importance of bare ice albedo and the systemic bias in GrIS surface mass loss that can result from neglecting spatially and temporally varying bare ice albedo, and the improved model physics that can be utilized to improve future projections of GrIS melt and sea-level rise.Bare ice albedo is kept constant in the visible (0.6) and NIR (0.4) bands of the default version of E3SM and our results indicate that such values are generally too high-with respect to the ones measured by MODIS-promoting an underestimation of surface melt.Other ESMs also use constant ice albedo values, such as CESMv2, which assumes 0.5 in the visible and 0.3 in the NIR (van Kampenhout et al., 2020).The CESMv2 bare ice albedo values are low biased compared to the average exposed ice area weighed albedo values presented here (0.56 and 0.33 respectively).As a result, the albedo and SMB biases presented in this study are not directly applicable to CESM.Future work would need to be done comparing the default CESMv2 bare ice scheme to our MODIS informed SNICAR-ADv4 technique.Future work should also compare our SNICAR-ADv4 enabled E3SM simulated SMB to measurements or regional models such as MAR or RACMO to quantify how strongly the SNICAR-ADv4 ice albedo scheme improves modeled SMB.However, previous work has shown that regional models simulate GrIS SMB more accurately than Earth system models as regional models have been calibrated to accurately simulate snow processes on the GrIS and have a much higher spatial resolution (Fettweis et al., 2020).This highlights the importance of the spatial resolution and sub-scale processes in global Earth System models.
The ability to use physically based radiative transfer calculations to simulate bare ice albedo on the GrIS is increasingly important, as exposure of bare ice on the GrIS is expanding both spatially and temporally (Ryan et al., 2019;Shimada et al., 2016;Tedesco et al., 2016).The improved model physics presented in this study allows ELM to more accurately simulate the fluctuating albedo on the GrIS, as opposed to simply prescribing MODIS bare ice albedo within the model.The implementation of SNICAR-ADv4 described here utilizes satellite retrievals of ice properties and therefore presents a realistic present-day climatology of the GrIS.Although the For the difference, the color scheme is the same as that in Figure 9c.
bare ice properties are limited to the MODIS era, the climatology of the bare ice properties can also be applied to both historical and future simulations.While the bare ice properties are inferred from 21st-century data, the combination of improved model physics and more realistic bare ice albedo is an improvement over the commonly used constant ice albedo schemes.This scheme allows for variation in ice albedo with varying insolation and atmospheric conditions and the quantification of sensitivities to ice physical properties and LACs.An avenue for future work is to utilize the MODIS inferred climatology to constrain the bare ice properties at higher elevations on the ice sheet, to be used as the snowline elevation increased in projections simulations.However, the use of prescribed ice properties in the model may inhibit important feedbacks, such as those involving temperature and SSA.The inability to simulate these feedbacks would likely impact climate simulations over longer timescales and future projections.Although, this could be remedied with an improved physical representation of ice properties within E3SM, as the SNICAR-ADv4 model can utilize evolving physical properties of ice that are prognosed in future versions of the model.Including evolving ice physics in ELM could be another avenue for future work.There is also work being done utilizing hyperspectral remote sensing techniques to better constrain snow and bare ice physical properties that can be utilized to better constrain LACs and snow and ice microphysical properties (Bohn et al., 2021;Vandecrux et al., 2022;S. Wang et al., 2018).
SNICAR-ADv4 enabled ELM can model the influence of bare ice exposure as it can now accurately simulate bare ice albedo while still capturing the effect of surface snow coverage and various LACs (Dang et al., 2019;Flanner et al., 2021;Whicker et al., 2022).While this work greatly improves our ability to more realistically simulate GrIS bare ice albedo, it does not provide information about the LAC content on the ice sheet.In the dark-ice region in the southwest ablation zone of the GrIS, the presence of LACs such as dust, black carbon, and darkly pigmented glacier algae heavily influence the albedo and contribute a significant amount toward total melt in the dark zone (Bohn et al., 2022;Chevrollier et al., 2022;Cook et al., 2020;Ryan et al., 2019;Tedesco et al., 2016;Tedstone et al., 2020).The current "equivalent black carbon" technique used to account for LACs allows SNICAR-ADv4 to accurately simulate bare ice albedo and surface melt, but it does not provide a robust analysis of the contributions from various LACs and their individual impacts on surface melt.

Figure 1 .
Figure 1.Schematic of SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) land ice column (right, utilized in this study) contrasted with the ELM default land ice treatment (left).SNICAR-ADv4 utilizes the land ice layers in ELM and prescribes the optically relevant physical properties of ice by inverting Moderate Resolution Imaging Spectroradiometer albedo data (MCD43A3 band 2 and band 8).The new bare ice albedo varies with the solar zenith angle, ice physical properties, and the presence of light absorbing constituents, in this case of this study that is limited to black carbon.

Figure 2 .
Figure 2. Albedo in the Moderate Resolution Imaging Spectroradiometer MCD43A3 band 2 region of the spectrum (0.841-0.876μm) as a function of ice specific surface area (SSA) and solar zenith angle.The relationships presented here are used to infer spatially and temporally varying bare ice SSA to prescribe in ELM.

Figure 3 .
Figure 3.The bare ice physical properties inferred from the inverse modeling technique, (a) specific surface area, (b) ice density, (c) air bubble effective radius.These spatio-temporally varying ice properties are utilized in the SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) enabled ELM simulations.

Figure 4 .
Figure 4. (a) Moderate Resolution Imaging Spectroradiometer (MODIS) bare ice near-infrared albedo and (b) offline SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) bare ice near infrared albedo on the ELM 0.5 × 0.5-degree grid.The offline SNICAR-ADv4 near-infrared albedo is generated in the offline lookup table simulations.

Figure 5 .
Figure 5. (a) Moderate Resolution Imaging Spectroradiometer (MODIS) bare ice visible albedo, (b) offline SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) bare ice visible albedo with no light absorbing constituents, (c) offline SNICAR-ADv4 bare ice visible albedo with equivalent black carbon concentrations determined by the offline lookup table and MODIS visible albedo, and (d) the inferred equivalent black carbon concentration.All plots contain the bare ice land area-weighted average in the bottom right corner.The offline SNICAR-ADv4 albedo values are generated in the offline lookup table simulations and the ELM area and bare ice fraction are used to weight the averages.

Figure 7 .
Figure 7.The ELM simulated average fraction of the grid cell that contains exposed bare ice in June, July and August.

Figure 8 .
Figure 8.A time series of the SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) Moderate Resolution Imaging Spectroradiometer-informed (line plot) and default ELM (dot plot) exposed bare ice weighted average (a) shortwave broadband albedo, (b) visible albedo, and (c) near-infrared albedo.The default ELM scheme overestimates the bare ice albedo in all cases except in June in the visible part of the spectrum.

Figure 9 .
Figure 9.The 2000-2021 annual average Greenland Ice Sheet wide climatic surface mass balance for (a) the default ELM simulation, (b) the SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) enabled simulation, (c) and the difference between the default and the SNICAR-ADv4 enabled simulations.In the surface mass balance (SMB) plots, red values indicate a mass loss while blue values indicate an increase in mass.In the difference plot, blue values indicate that the default scheme has higher SMB values, and the red values indicate that SNICAR-ADv4 scheme has larger SMB values.

Figure 10 .
Figure 10.The summertime average Greenland Ice Sheet wide climatic surface mass balance for (a) the default ELM simulation, (b) the SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4 (SNICAR-ADv4) enabled simulation, (c) and the difference between the default and the SNICAR-ADv4 enabled simulations.For the difference, the color scheme is the same as that in Figure 9c.

Table 1
The Exposed Bare Ice Area-Weighted Average June, July and August Albedo in the Shortwave, Visible, and Near-Infrared Bands When the Moderate Resolution Imaging Spectroradiometer (MODIS) values are weighted by the actual exposed bare ice regions, they are lower than that prescribed in the default ELM simulations.The SNICAR-ADv4 enabled simulations show good agreement with the MODIS measurements.Abbreviation: SNICAR, SNow, ICe and Aerosol Radiation model Adding-Doubling Version 4.