Evaluation of a new snow albedo scheme for the Greenland ice sheet in the regional climate model RACMO2

. Snow and ice albedo schemes in present day climate models often lack a sophisticated radiation penetration scheme and do not explicitly include spectral albedo variations. In this study, we evaluate a new snow albedo scheme in the regional climate model RACMO2 for the Greenland ice sheet, version 2.3p3, that includes these processes. The new albedo scheme uses the two-stream radiative transfer in snow model TARTES and the spectral-to-narrowband albedo module SNOWBAL, version 1.2. Additionally, the bare ice albedo parameterization has been updated. The snow and ice broadband and narrow- 5 band albedo output of the updated version of RACMO2 is evaluated using PROMICE and K-transect in-situ data and MODIS remote-sensing observations. Generally, the modeled narrowband and broadband albedo is in very good agreement with satellite observations, leading to a negligible domain-averaged broadband albedo bias smaller than 0.001 for the interior. Some discrepancies are, however, observed close to the ice margin. Compared to the previous model version, RACMO2.3p2, the broadband albedo is considerably higher in the bare ice zone during the ablation season, as atmospheric conditions now alter 10 the bare ice broadband albedo. For most other regions, however, the updated broadband albedo is lower due to spectral effects, radiation penetration or enhanced snow metamorphism.


Introduction
The absorption of shortwave radiation is an important component of the surface energy balance of snow-covered surfaces (Van den Broeke et al., 2005;He et al., 2018b;Warren, 2019). A drop in the surface reflectivity for solar radiation, i.e., albedo, 15 leads to more absorbed energy in the snowpack, which in turn leads to higher snow temperatures or melt. This melt-albedo feedback is initiated once snow starts to melt and the snow structure is altered, lowering the albedo (Van As et al., 2013;Jakobs et al., 2019). It is therefore imperative for regional and global climate models (RCMs and GCMs, respectively) to capture snow albedo correctly in order to reproduce the current climate and to make future climate projections for snow covered glaciated regions such as the Greenland ice sheet (GrIS). 20 The albedo of snow and ice is highly spectrally dependent and also depends on various other quantities. For clean snow, the spectral albedo, i.e. the albedo as a function of wavelength, is almost one for near-ultraviolet (near-UV, 300-400 nm) and visible light (400-750 nm), but drops for near-infrared (near-IR, 750-1400 nm) and is low and fluctuating for infrared (IR) radiation (Fig. 1d, Warren and Wiscombe (1980); Warren et al. (2006); Gardner and Sharp (2010); Dang et al. (2015); Picard et al. of the parameterization of Gardner and Sharp (2010) to derive a broadband albedo. Therefore, RACMO2 until now did not include explicit spectral albedo or spectral irradiance effects, nor an adequate radiation penetration scheme. Introducing a new snow albedo parameterization that includes these processes is therefore timely.
RACMO2.3p3 uses a new snow albedo parameterization using the Two-streAm Radiative TransfEr in Snow model (TARTES, Libois et al., 2013) coupled with the Spectral-to-NarrOWBand ALbedo (SNOWBAL) module version 1.2 (Van Dalum et al., 2019). This set-up provides appropriate narrowband albedos for all 14 shortwave spectral bands utilized in RACMO2. TARTES also considers radiation penetration for its surface albedo calculations and provides estimates of energy absorption in the snowpack. Additionally, the new snow albedo parameterization is used to update the :::::: develop : a :::: new : ice albedo scheme.
Here, we present and evaluate the broadband and narrowband albedo modeled by RACMO2.3p3 for the GrIS, and compare it with remote sensing data, in-situ observations and the broadband albedo modeled by the previous iteration of RACMO2, ver-70 sion 2.3p2. The remainder of this manuscript is made up of six sections. Section 2 summarizes the changes made in RACMO2 and introduces the remote sensing and in-situ observational data sets. Section 3 and Sect. 4 evaluate the new RACMO2 version with these data sets. Comparisons between the albedo modeled by RACMO2 version 2.3p3 and 2.3p2 are shown in Sect. 5.
Finally, the sensitivity to the chosen impurity concentration of snow is analyzed in Sect. 6, and results are discussed and conclusions are drawn in Sect. 7. The impact of the model improvements on the climate, surface mass balance and surface energy 75 balance of the GrIS ice sheet will be discussed in a forthcoming publication.
2 Model and observational data sets 2.1 Regional climate model version of RACMO2, version 2.3p2, from now on abbreviated to Rp2, is adapted for glaciated regions :::: tiles by using a multilayer snowpack that interacts with the atmosphere and involves processes within the snow column, such as melt and refreezing. Rp2 is introduced in more detail in Noël et al. (2018). At the lateral boundaries, RACMO2 is forced with ERA-Interim data (Dee et al., 2011). In the new RACMO2 version, 2.3p3, from now on abbreviated to Rp3, two components have been adjusted, the 85 multilayer firn module and the snow and ice albedo parameterizations for glaciated regions.
Firstly, Rp2 uses :::: used a prognostic fresh snow layer, which is effectively a sublayer of the uppermost snow model layer. In Rp3 this fresh snow layer is removed; instead, the uppermost layers are allowed to be very thin, i.e., in the order of millimeters, thus containing fresh snow only. For heat diffusion calculations, these thin layers are treated as a single layer to maintain numerical stability. If melt or refreezing occurs in one of these layers, their individual temperature is estimated obeying the 95 temperature gradient and conserving their combined heat content.
Secondly, in Rp2 layers below a threshold thickness merge :::::: merged : with the first layer below. In Rp3, a layer merges with the adjacent layer having most similar density and grain size. Furthermore, undesired numerical diffusion is avoided by implementing mass redistribution if a thin layer merges with a thick layer. A layer containing glacial ice is not allowed to merge with layers that are formed locally, i.e., by snow deposition on this grid point. This allows for the formation and preservation 100 of layers with ice lenses.
Melting of the uppermost layer always leads to thinning, regardless of its density.

Snow albedo
Rp2 uses :::: used a plane-parallel broadband snow albedo scheme based on Gardner and Sharp (2010), that depends :::::::: depended indirectly on wavelength in the form of tuning parameters, and is limited to two snow layers. This albedo scheme parameterizes 115 ::::::::::: parameterized : albedo variations due to a changing SZA, grain radius, cloud cover, impurities, and an altitude-dependent atmospheric optical thickness, the latter for clear-sky conditions (Kuipers Munneke et al., 2011). In RACMO2, the first two snow layers are often very thin, i.e. a few millimeters for fresh snow and up to a few centimeters for older snow, thus effectively neglecting almost all radiation penetration.

140
For the other bands, three narrowband albedos are determined, i.e., for direct and diffuse radiation for clear-sky conditions, and for diffuse radiation for cloudy conditions. Clear-sky and total-sky narrowband and broadband albedos are then determined using the modeled radiative fluxes. Note that clear-sky and total-sky albedo are identical if no clouds are present. Finally, for evaluation with the seven MODIS narrowband albedos, clear-sky diffuse radiation albedos are also explicitly derived for these bands.
Next, the MODIS bare ice albedo range is converted to an impurity concentration. Using the fitted SSA, ice density and a reference SZA of 60 • , which is the largest angle for which the observations of MODIS for the GrIS are still somewhat reliable (Wang and Zender, 2010), together with RACMO2 narrowband irriadiance profiles for such conditions, a broadband albedo 175 can be calculated for a range of impurities such that the MODIS albedo range is covered. The resulting impurities vary :::: soot ::::::::::: concentration ::::: varies : between 69 ng g −1 for an albedo of 0.55, to 2445 ng g −1 for an albedo of 0.3 (Fig. 1c), and are saved in a lookup table. This lookup table is used to convert the MODIS albedo field, after resampling to the 11-km grid of RACMO2, to an impurity field feasible for TARTES to use in RACMO2 (Fig. 1b). Adding soot alters, as expected (Doherty et al., 2010;Gardner and Sharp, 2010;He and Flanner, 2020), only the spectral curve in the near-UV and visible part of the spectrum (Fig.   180 1d). Note that the broadband albedo can still reach values beyond the range indicated in Fig. 1c, depending on atmospheric conditions and SZA.

In-situ measurements
In-situ observations provide insight in the performance of RACMO2 for total-sky and cloudy conditions, unavailable from remote sensing observations. Therefore, we evaluate the albedo product with the Automatic Weather Station (AWS) data along PROMICE AWS are mostly located in the ablation zone. Only AWS sites that cover at least partially the model period and for which an appropriate grid point can be selected in RACMO2, are selected. Figure 2c shows the location of the PROMICE

285
Most pronounced biases can be observed for grid points close to the margin in the southeast (region D), but extending beyond this region. The albedo difference is typically around 0.03, but can be as high as 0.15. The ablation zone in D is too narrow (up to 10 km) to be adequately resolved at the used resolution of 11 km (Noël et al., 2015). Furthermore, the southeast is characterized by heavy snowfall (Ohmura and Reeh, 1991;Mernild et al., 2015), and in combination with the low model resolution, snow persists at the surface throughout the ablation season in Rp3. Still, the grain radius increases and the albedo 290 drops somewhat during the ablation season, but not as fast as the albedo decline in MODIS. Note that the uncertainty of MODIS is also considerable, as this is mountainous terrain (see Sect. 2.3). respectively, instead of bare ice that MODIS observes. These data points originate from the southeast (area D of Fig. 2c, see the previous discussion). The bias for the Summit region (Fig. 3b) is similar to both the bias of area A (Fig. 2c) and the bias of the MODIS WSA product (Stroeve et al., 2013) :::: CSD :::::: albedo ::::::: product ::::::::::::::::: . Furthermore, the albedo variability 300 is low for this region, as melt does normally not occur here (except for July 2012 (Nghiem et al., 2012;Bennartz et al., 2013)) and metamorphosis :::::::::::: metamorphism : is slow in this cold climate. Figure 3c represents an area in southern Greenland without bare ice, which is characterized by albedo decrease due to rapid snow metamorphism (Lyapustin et al., 2009). Rp3 performs well for this region. The region shown in Fig. 3d corresponds roughly with the ablation zone of the K-transect. The bias is very indicating that it is not the bare ice albedo scheme that causes the discrepancies.

325
Along the K-transect, albedo data are available at S5, S6, S9 and S10 (Fig. 2c) for 15:00 UTC and are shown in Fig. 5 for 2012 and compared with Rp3 clear-sky and total-sky conditions, Rp2 total-sky, and MODIS WSA :::: CSD. S5 is too close to the ice margin for MODIS to produce a reliable albedo at this resolution, and is not included. As the area around this station is characterized by very rough terrain, RACMO2 at this low resolution has trouble reproducing the in-situ albedo.
For S6, Rp3 albedo corresponds well with the measurements, especially when snow returns after the melt season. At the 330 onset of the accumulation season, Rp2 performs :::::::: performes : considerably worse, as a lack of radiation penetration causes a too rapid albedo increase. As the accumulation season progresses and the snow layers become thicker, the albedo difference between Rp3 and Rp2 diminishes. MODIS WSA :::: CSD :::::: albedo fits well with Rp3 clear-sky albedo for bare ice conditions. The resolution of RACMO2 is not sufficient to capture all spatial variations, i.e., the AWS data might not be representative for the full grid point of RACMO2 in which S6 is located.

335
For S9 (Fig. 5c), both Rp3 and Rp2 total-sky albedo show relatively large deviations with the in-situ measurements during the accumulation season, but the difference of Rp3 clear-sky albedo with MODIS WSA :::: CSD :::::: albedo is much smaller. Site S9 is characterized by spatial inhomogeneity, although to a lesser degree than S5 and S6. During the start of the ablation season, superimposed ice persists at the surface, delaying the albedo drop to bare ice values by a few days. Just like at S6, the Rp3 albedo fits better than Rp2 with both in-situ and MODIS observations at the onset of snowfall after summer. 340 S10 observations (Fig. 5d), which are representative for the lower accumulation zone of Greenland, correspond well with both Rp3 and Rp2. In addition, Rp3 fits with MODIS WSA :::: CSD :::::: albedo : before the melt season, and performs reasonably well during and after the melt season. 2012 is characterized by a long and intense melt season, explaining the albedo decrease observed in MODIS WSA :: the ::::::: MODIS ::::: CSD :::::: albedo ::::::: product. The melt season albedo decrease in Rp3, however, is slightly delayed.

14
In this section, we compare the Rp3 albedo product with Rp2 and highlight the differences. Moreover, we investigate the impact that clouds have on the albedo and investigate the seasonal differences. Figure 6 shows the 16-days running mean total-sky albedo for 15:00 UTC between 2006 and 2015 for Rp2 and the albedo 350 difference between Rp3 and Rp2. For most of the ice sheet, a small negative difference is observed, i.e., the albedo of Rp3 is slightly lower than Rp2. In areas E in the north-west and F in the north-east (Fig. 6b), the albedo of Rp3 is lower than Rp2. Both regions are characterized by limited snowfall, resulting in a shallow snowpack on top of bare ice during the accumulation season for several months (Fig. 7). For September and October, while the sun is still above the horizon, snow layers are thin enough 355 for solar radiation to reach the bare ice layersunderneath, :: At ::: the ::: end ::: of ::: the ::::::: ablation :::::: season ::::: when :::: new ::::: snow ::::: layers :::: start ::: to ::::
We also observe a small albedo difference (smaller than 0.01) in the high accumulation region in southeast Greenland (area D of figure 2c). For some grid points around the margins, the albedo of Rp3 is considerably higher than Rp2. This difference is more pronounced in areas with exposed bare ice during the ablation season, and is limited to a single grid point bordering nonglaciated tiles. These differences can be traced back to uncertainties in the bare ice albedo field of Rp2, where grid points close 365 to the margin are :::: were contaminated with tundra albedo, and an albedo difference with the new model is therefore expected.  2c. In general, Rp3 albedo correlates well with the albedo product of Rp2 (Fig. 8a), as the bulk of occurrences are close to the 1-to-1 line (black line). Most of the GrIS is covered in snow, for which the albedo is high (larger than 0.75) and the albedo difference is small (number 1 in Fig. 8a). Snow metamorphism is slightly stronger in Rp3, leading to lower albedos (about 370 0.75) compared to Rp2 (about 0.8, 2). Both 1 and 2 typically occur in the interior (Fig. 8b).
Larger albedo differences occur when firn or ice are close to the surface as radiation penetration lowers Rp3 albedo (3). A snow profile similar to Fig. 7 in early June would result in such an albedo difference. For the 0.6 albedo bin of Rp2 (4), the albedo of Rp3 is generally 0.05 higher, illustrating that Rp3 seems to have a slower firn-ice transition. In south Greenland (Fig.   8c), processes 1, 2, 3 and 4 are all relevant.

375
The snow and ice albedo merge well, but the Rp3 albedo is often higher than Rp2 for the ablation zone, i.e., for albedos smaller than 0.6. Still, highest occurrence density is found near the 1-to-1 line (5). As the new bare ice albedo parameterization allows variations due to atmospheric conditions, as is described in Sect. 2.1.3, some deviations are expected. That is to say, higher albedos with respect to Rp2 are expected for bare ice, as atmospheric variations usually increase the albedo, e.g., clouds.
Additionally, edge errors in Rp2 cause :::::: caused a considerable albedo difference, typically when Rp2 albedo is :::: were : smaller than 0.35 (6). Processes 1 to 6 are all applicable to the region around the K-transect (Fig. 8d), causing a spread for this region that is considerably larger than for the regions in Fig. 8b and Fig. 8c.

390
To conclude, the albedo product of Rp3 is often similar to Rp2, but some distinct differences are still notable ( Fig. 6 and Fig.   8). All these differences can be well understood in terms of physical processes.
As clouds absorb mainly in the (infra)red part of the spectrum, a blue shift occurs, which consequently increases the :::::: surface 395 albedo (Dang et al., 2015). The cloud dependence of the ::::: surface : albedo parameterization in Rp2 is ::: was : limited to the cloud optical thickness, neglecting water vapour while no distinction is made between liquid or ice water clouds. Figure 9 shows for a grid point in south-central Greenland (star in Fig. 6b) the total-sky ::::: surface : albedo, transmissivity for shortwave radiation in the atmosphere, i.e., the ratio between the top-of-atmosphere (TOA) and surface downwelling radiation, and the fraction of TOA shortwave radiation absorbed in snow as a function of vertically integrated cloud content (VICC), i.e., total liquid water 400 and ice in the atmosphere above a point, for Rp3 and Rp2.
For clear-sky conditions (VICC smaller than 0.05 kg m −2 ), the :::::: surface albedo in Rp3 is slightly lower than Rp2 (Fig. 9), which is in agreement with previous results. The :::::: surface : albedo in Rp3, however, increases more rapidly with VICC than in Rp2, leading to higher albedos for large VICC. Furthermore, the transmissivity decreases slower with VICC than in Rp2 and less radiation is absorbed in snow. For example, for a VICC of 0.3 kg m −2 , the :::::: surface broadband albedo, transmissivity 405 and fraction of absorbed energy changes by 0.02, 0.04 and -0.01, respectively. The latter implies that the net SW absorption decreases by approximately 25%. These differences show that as a cloud thickens and the :::::: surface albedo increases in Rp3, more reflected radiation will interact with clouds, eventually raising the transmissivity. Consequently, as these shorter wavelengths now scatter more often and are less likely to be absorbed in the snow, a white-out effect occurs, which is not captured in Rp2.
In a Taylor diagram, the azimuthal position illustrates the correlation coefficient, the radial distance to the origin the standard deviation and quarter circles with its origin at the 1.0 standard deviation the misrepresented variability. A data set matches the observations perfectly if it is located at the star in Fig. 12a, i.e., with a correlation coefficient and standard deviation of 1 and a misrepresented variability of 0. Data sets located away from the star but on the dashed line have the same variance as the 445 observations, but do not correlate perfectly, leading to a higher misrepresented variability. Data sets close to the origin and close to the y axis, for example, are characterized by an underestimation of the standard deviation and a low correlation coefficient, while data sets beyond the dashed line and close to the x axis overestimate the standard deviation, but have a high correlation coefficient. Similarly, the misrepresented variability is illustrated on the y axis, bias on the x axis and the RMSE on semi circles in Fig. 12b, with data sets performing better close to the origin.

450
Model performance varies for each observational site, but the sensitivity to small impurity concentrations is generally low.
High concentrations of impurities such as 50 ng g −1 , on the other hand, alter the albedo considerably, reducing the quality of Rp3 for almost all stations. For NUK-U and QAS-U, which are both located in south Greenland and within 50 km of the ice margin, RACMO2 correlates relatively well, but severely underestimates variability and shows a large bias. As the ice margins are characterized by a high soot concentration (Doherty et al., 2010), a higher modeled soot concentrations consequently 455 performs somewhat better for these locations. Snow cover in RACMO2 is also too homogeneous and similar processes that we discussed for area B and D of Fig. 2c occur. Additionally, RACMO2 is known to overestimate snowfall for QAS-U, inhibiting bare ice to surface (Noël et al., 2018). For S10 and KAN-M, which are most representative for the interior of the ice sheet, the differences between Rp3 with impurity concentration of 5 and 10 ng g −1 and to a lesser degree also with 0 ng g −1 and Rp2, are small and in good agreement with observations. To summarize, as the sensitivity of Rp3 to small impurity quantities is low 460 and the snow albedo is generally in good agreement with observations, an ice-sheet wide impurity concentration of 5 ng g −1 is a safe choice and for the interior is in good agreement with observations (Fig. 11b).

Summary and conclusions
We evaluated the new spectrally-dependent snow and ice albedo parameterization in RACMO2, based on TARTES and coupled by SNOWBAL, for the Greenland ice sheet. The albedo correlates well with the MODIS MCD43A3 white-sky ::::::: clear-sky :::::: diffuse 465 albedo product for both broadband and its seven spectral bands, and performs especially well in the interior. Some discrepancies around the margins are observed, which can be partly ascribed to resolution problems and excessive modeled snowfall, but also to uncertainty in the MODIS product. Around the K-transect, for which many observations are available, the snow and ice albedo in RACMO2 show acceptable small deviations with in-situ and MODIS observations.
With respect to the previous albedo parameterization of RACMO2, slightly lower broadband albedos are modeled. Although 470 large broadband albedo differences at the margins are due to an error in the old version, most changes can be ascribed to improved physics. Radiation penetration, subsurface heating, the inclusion of narrowband albedo and spectral shifts due to solar zenith angle, water vapour and both ice and water clouds are now all incorporated.
There is, however, still room for improvement. The soot concentration for snow is fixed in RAMCO2, while it can change considerably over space and time (Chylek et al., 1992;Doherty et al., 2010;Van Angelen et al., 2012;Dang et al., 2015).

475
Although RACMO2 shows a low sensitivity to small impurity concentrations, a prognostic soot model for snow prescribing a dynamic one-dimensional soot concentration profile is still preferable. Furthermore, no other impurity types are included.
We have also improved the bare ice albedo field by coupling it with TARTES and defining a fictitious grain size, which allows the broadband bare ice albedo to vary with atmospheric conditions. Nevertheless, a proper bare ice albedo module would still be preferable, as for example, a process like darkening due to biological activity is not incorporated, which may become the Taylor diagram, which shows the misrepresented variability on the y axis, bias on the x axis and RMSE on semi circles from the origin.
Evaluation of the narrowband albedo of RACMO2 remains limited, as most of the spectral bands of MODIS are located at the edge of rather large bands of RACMO2 and hence cannot be compared directly. To solve this, we have run TARTES specifically for diffuse radiation for the seven bands of MODIS. For these bands, differences between modeled albedo and 485 observations are low. In the bare ice zone, band 1, 3 and 4, which are within the visible light part of the spectrum, show a larger narrowband albedo bias than the other bands. Larger spatial variations are observed for bands 2 and 5, which are characterized by a strong sub-band spectral albedo gradient.
To conclude, the new snow and ice albedo scheme of RACMO2 performs very well compared to remote sensing and insitu observations for the Greenland ice sheet. Differences with the previous :::::::: RACMO2 : version are generally small, but where differences are observed, the new processes lead to improved broadband albedo estimates. The improvement of the albedo scheme of RACMO2 will enhance its ability to make future climate projections. In a forthcoming publication, we assess the impact of the new snow and ice albedo scheme on the surface mass balance, surface energy balance and subsurface energy absorption of the Greenland ice sheet.
Author contributions. CTvD, WJvdB and MRvdB started and decided the scope of this study and analyzed the results. CTvD implemented 495 the new snow albedo scheme in the model, performed the simulations and led the writing of the manuscript. SL processed and provided remote sensing data. All authors contributed to discussions on the manuscript.
We would like to thank the reviewers for their constructive comments that have improved the accuracy of the evaluation and the clarity of the paper. In black the comment, in orange the response, in blue the changes that we would implement in the manuscript.

Review #1 by Stephen Warren
(1) Lines 28-29. "With coarser grains, the likelihood for light to reflect off a grain's surface out of the snowpack reduces, lowering the albedo." This is not true; reflection off a grain's surface is independent of grain size and is anyway a minor contributor to the albedo; successive refraction instead dominates (Bohren and Barkstrom, 1974). The correct explanation for dependence of albedo on grain size is given by Warren (2019, cited in the paper): "In coarse-grained snow, a photon travels a longer distance through ice between opportunities for scattering than in fine-grained snow, so it is more likely to be absorbed, and therefore a snowpack of coarse grains has lower albedo." As our description was indeed inaccurate, we changed the following: Page 2: With coarser grains, light has to travel longer through ice before it has the opportunity to reflect off a grain's surface out of the snowpack than for fine-grained snow, hence lowering the albedo (Wiscombe and Warren, 1980;Gardner and Sharp, 2010;Picard et al., 2012, Warren, 2019. Fresh snow with a small grain radius… (2) Lines 140-142. "We assume that clean bubble-free ice has an albedo of approximately 0.6 (Reijmer et al., 2001) and that the bare ice albedo is subsequently lowered by . . . bubbles . . . ." This statement is wrong in three ways. First, Reijmer et al. did not say that bubble-free ice has albedo 0.6. Second, a thick layer of ice, if it is really bubble-free, will have very low albedo, about 0.07. Bubblefree ice does sometimes occur on frozen lakes, where it is appropriately called "black ice". Third, the albedo of ice is raised by bubbles, not lowered. The bubble surfaces are what are responsible for the model's assumed specific surface area of 0.788 m 2 /kg. Without bubbles (or cracks), the SSA of thick ice would be zero, and the albedo would be only the Fresnel reflection at the top surface (which is 0.07 for diffuse incidence). Figure 17b of Dadic et al. (2013, cited in the paper) shows the broadband albedo versus SSA for ice, firn, and snow, all on the same plot. You are right. Therefore we changed the following: Page 5: Firstly, we assume that clean blue ice has an albedo of approximately 0.6 (Reijmer et al, 2001;Dadic et al, 2013). Blue ice is typically found in areas with a very smooth surface and high sublimation rates, but no melt, and has a high bubble content, leading to a relatively high albedo. The bare ice albedo is subsequently lowered by standing water, surface roughness and impurities. Furthermore, we assume that MODIS… (3) Superimposed ice is mentioned in several places (e.g. lines 285-286), but is not discussed adequately. The reader wants to know the definition of superimposed ice, and why its albedo is higher than that of bare glacier ice. Glacier ice is formed by compression of snow under pressure, whereas superimposed ice is formed by refreezing of meltwater. The paper describes superimposed ice as having a granular structure with grain radius 0.7-1.0 mm, which could result in a higher SSA than in glacier ice. But there are reports of such a surface granular layer developing also in glacier ice exposed to sunlight, as in Figure 1 of Mueller and Keeler (1969). We agree that we have not clarified the definition of superimposed ice clearly. Exactly as you have stated, we mean with superimposed ice layers that have refrozen meltwater in it, which is different from bare ice that is formed by compaction. More importantly for us, however, is the different way that superimposed ice is treated compared to bare ice in the model. Superimposed ice does not use the bare ice albedo scheme, as is described in Sect. 2.1.3, but uses the snow albedo scheme and the changes described in Sect. 2.1.4. We think that changing the following in Sect. 2.1.4 is sufficient for the reader to know what we mean with superimposed ice for the rest of the paper. Page 6: In Rp3, superimposed ice is treated differently than glacial ice. Superimposed ice forms in snow layers by refreezing of melt water, while glacial ice forms by compaction of snow. As superimposed ice has a granular structure , it has to be treated differently than bare ice. Due to the granular structure, it is desirable to use the snow albedo scheme of Sect. 2.1.2 over the bare ice albedo scheme of Sect. 2.1.3, as a fixed rather large grain radius is used for bare ice. However, without additional corrections, the typical grain size of model layers with superimposed ice in Rp3 is 0.7 to 1.0 mm, leading to unrealistically high albedos of approximately 0.7. The typical albedo of superimposed ice is 0.65, as is measured by  at S9 of the K-transect. In order to improve this, a minimum grain radius for superimposed ice is imposed, which increases linearly from 0.720 mm for a density of 750 kg m 3 , to the bare ice value of 4.152 mm for a density of 917 kg m 3 . This superimposed ice layer uses the same impurity concentration as a snow layer. This correction leads to realistic albedos for superimposed ice and exposed ice lenses.
(4) The effect of a thin snow layer on top of ice is belittled in the Rp3 model. On Line 329, Location 8, with an Rp3 albedo of only 0.45, is said to be thin snow on bare ice. But in a comparable situation, Figure 2 of Brandt et al. (2005) showed that when thick bare sea ice of albedo 0.49 was covered with 5-10 mm of snow, its albedo rose to 0.81. So why was the Rp3 albedo so low? You might check to see whether the snow was melting. After analysis of these cases, we agree that our conclusion was not correct. Location 8 represents grid points where Rp3 models either bare ice or melting snow or firn, while Rp2 models a fresh snow top layer. So it does not show a strong direct radiation penetration effect on albedo. We have changed the following: Page 13: At the end of the ablation season when new snow layers start to form, melt can still occur, forming melt water within those layers. This is more common in Rp3, as internal heating increases the subsurface temperature and therefore more melt occurs, changing the snow structure and lowering the albedo in Rp3 with respect to Rp2. Furthermore, if melt is strong enough, it can remove or thin the fresh snow top layer, exposing bare ice underneath. As the accumulation season… Page 15: … and to a lesser extent area E and some grid points close to the ice-sheet margin in northern and eastern Greenland. These occurrences represent cases that Rp2 models a fresh snow top layer without melt while Rp3 models either bare ice or a melting snow or firn layer.
Minor comments Line 32. In spite of the importance of snow grain shape for albedo, a model of spherical grains still accurately reproduces the measured spectral albedo, by adjusting the grain radius. The reason this procedure works was explained by . The reviewer is right and we have added the following in the manuscript. Page 2: …out of the snowpack (Libois et al., 2013;He et al., 2018a), but  show that a model with spherical grains can still accurately reproduce the measured spectral albedo by adjusting the grain radius. To summarize, it is thus essential… Line 38. "Rayleigh scattering [by the atmosphere] is more effective for larger wavelengths…". Change "larger" to "shorter". Rayleigh scattering is inversely proportional to the fourth power of the wavelength, so blue is scattered more than red. Changed as requested Page 2: …as Rayleigh scattering by the atmosphere is more effective for shorter wavelengths… Line 42. Change "incoming radiation" to "incoming direct-beam radiation". Done Page 2: …of the incoming direct-beam radiation… Line 54. "version 2.3p3". I am guessing that the "p" means "polar". Is that correct? Yes, "p" means polar. To clarify this, we have added the following: Page 2: The polar (p) version of RACMO2 is a model developed..
Line 64. "the new snow albedo parameterization is used to update the ice albedo scheme". Does this mean that snow and ice both use the same parameterization, just with different coefficients? Snow and ice now both use the TARTES/SNOWBAL albedo parameterization, but the grain radius/SSA and impurity content for bare ice is different. Page 3: Additionally, the new snow albedo parameterization is used to develop a new ice albedo scheme.
Lines 77-78. "The polar version of RACMO2 . . . is adapted for glaciated regions by using a multilayer snowpack . . ." Is this scheme used also for seasonal snow in midlatitudes, or only for polar latitudes? If the latter, what is the latitude domain to which the model is applied? Multilayer snowpack is not used for seasonal snow in midlatitudes, it is only used for glaciated tiles, which is for our domain the Greenland ice sheet and its surrounding glaciers, if those glaciers are large enough to fill up a grid point. Page 3: … adapted for glaciated tiles by using … Lines 97-98. "For ice layers . . . only a density reduction takes place and no thinning occurs." This is backwards. When ice melts, it thins, but the density of the remaining ice is unchanged. Paper Although you are right that the density of the remaining ice will remain, the density of the model layer will go down as pore space is created. Moreover, if melting ice reaches the surface, it will thin the layer. To clarify this, we changed the following: Page 4: Thirdly, internal energy absorption heats subsurface snow layers and can induce melt. In Rp3, melt will only thin a subsurface snow layer, i.e., a layer with a density below 700 kg m -3 and not change its density. For ice layers, i.e., with a layer density larger than 830 kg m -3 , melt creates pore space, reducing the layer density and no thinning occurs.
Line 104 (and elsewhere). "Rp2 uses . . ." To distinguish Rp2 and Rp3, I suggest using past-tense for Rp2 and present-tense for Rp3. As requested, we changed the tense to the past on relevant places when Rp2 is used.
Line 114. "grain shape . . . are all explicitly resolved". How is grain shape described in the model?
In the model, all grains are assumed to be spherical. Page 4: …are all explicitly resolved. In this study, all grains are spherically shaped. TARTES is able to calculate… Lines 145-149.  was the first one to show that bubbly ice could be modeled as verycoarse-grained snow, as done here. Thank you for this suggestion. We have added a reference to the work of  in the manuscript.
Page 5: …and thus can be used to indicate clean, bare ice, which is similar to the findings of .
Line 161. "unrealistically high albedos of 0.7". A reference is needed here for measured albedos of superimposed ice, to show that an albedo of 0.7 is too high. As is described in comment (3) of this review, we have added the following: Page 6: The typical albedo of superimposed ice is 0.65, as is measured by    Doherty et al. (2010) found an average of only 3 ng/g soot, but as the snow melts the soot accumulates at the surface, attaining concentrations in the top centimeter of 20-100 ng/g at Dye-2 ( Figure 8a of Doherty et al. 2013).
We have specified where this 3 ng/g is a representative number. We chose to have the impurity content right for non-melting conditions, in order to get the melt-onset timing right. Page 7: In the interior and when no melt occurs, the soot concentration is approximately 3 ng g -1 … Lines 178-179. "the bare ice albedo value is replaced with the bare ice soot concentration . . ." Rewrite this sentence. It does not make sense to replace an albedo value with a soot concentration; they have different units.
We have adjusted the formulation as the initial formulation resembled the numerical implementation of boundary files, while for the paper this is indeed only confusing. Page 7: …initialization day, i.e., 1 September 2000, but the fresh snow sub-layer data is omitted. Page 7: Furthermore, Rp3 is initialized with a soot concentration that can be used to calculate the bare ice albedo (Sect. 2.1.3).
Lines 233-234. "The region indicated by C . . . albedo difference during winter . .." Region C is at latitude 72.5 degrees, so on the March equinox at noon the SZA is 72.5 degrees, and throughout the winter the SZA is greater than this. Since you are limiting your comparisons to SZA<55 degrees, in winter the Sun is too low for any comparisons. You are right, we made a mistake with the analysis for region C, as the MODIS product cannot be used in winter for this location. The albedo for a grid point as a function of time is shown in the added figure (left) for a grid point within area C for 2012. Clearly, Rp3 underestimates the clear-sky albedo here. This difference is likely caused by the small amount of snow present in Rp3, which quickly thins and vanishes during summer, reducing the albedo (right figure). The region indicated by (C) shows an underestimation of Rp3 albedo compared to MODIS. Rp3 typically only models a few decimeters of snow during winter that melt quickly during summer, thinning and eventually removing these layers. Throughout the summer, Rp3 likely underestimates the snow thickness, leading to underestimated albedo.
Lines 252-253. "metamorphosis . . . metamorphism". Use consistent terminology. LaChapelle favored "metamorphism". Changed as requested Lines 269-270. "large quantities of soot . . . not enough to lower the albedo for bare ice to MODIS values". Bare ice on Greenland contains not just soot but also dust, cryoconite, and algae, all of which reduce the albedo.  , Table 1) reported ice albedos as low as 0.2. We have adjusted the text to include this notion. Page 12: It is, however, expected that Rp3 needs large amounts of soot to lower the bare ice albedo, as no dust, cryoconite and algae are modeled, all of which lower the albedo .
Line 335 Section 5.2. In this section, clarify that by "albedo" you mean surface albedo, not TOA albedo. As requested, we changed the term 'albedo' to 'surface albedo' in this section.
Line 366. "[In autumn] the SZA is too small to make a significant difference". The Sun is nearly as low in SON as in DJF; they are the two "low-sun" seasons. So what you said about the effect of SZA in DJF should also apply to SON. You are right, and we change it accordingly. Page 18: Autumn (September, October and November (SON) , Fig 10d) presents small deviations for most of the ice sheet, as fresh snow layers accumulate and the SZA increases, slowly transitioning to winter conditions and the processes described for DJF become increasingly important.
Line 417 states that the TARTES snow model is based on Rayleigh scattering. Rayleigh scattering applies to gases, and to particles much smaller than the wavelength. Maybe TARTES uses Rayleigh theory for submicron soot particles, but surely not for snow grains.
TARTES uses a geometric-optics method, so not Rayleigh scattering, so you are right. We have changed the following: Page 4: …asymptotic analytical radiative transfer theory (Kokhanovsky, 2004;He and Flanner, 2020) using the geometric-optics method… Page 5: Despite not using Mie-scattering theory, the spectral curve in TARTES for this SSA value resembles the expected curve for bare ice… Page 19: TARTES is based on simple geometric-optics theory, while an ice albedo module with Mie-scattering theory is… Figure 1. Much of the coastal region, but not all, is blank. What are the blank areas; are these ice-free and snow-free? That is correct, the blank areas are not glaciated. For the changes made in the text: see next comment Figure 1 caption. Line 1. "16-day". What months are included? Lines 1-2. "The albedo is limited between 0.30 and 0.55". I don't understand this; you know that the accumulation area has albedo >0.55. Line 3. "impurities". What kind of impurity is assumed, and what is its mass-absorption crosssection? As requested, we added some missing information and changed the caption of Figure 1 to: Page 6: Figure 1 (a) Lowest five percent of the MODIS MCD43A3v5 1 km 16-day clear-sky diffuse albedo field for glaciated areas for the period 2000-2015. As this albedo field is used to determine a bare ice albedo field, it is limited between 0.30 for dark ice in the ablation zone, and 0.55 in the accumulation zone under perennial snow for consistency with RACMO2 (Noël et al., 2018). (b) Bare ice impurity field that is implemented in RACMO2.3p3 for glaciated grid points. Here, all impurities are soot. (c) In blue, fitting the specific surface area (SSA) to clean blue ice albedo, which is assumed to be 0.6. The fitted SSA equals 0.788 m 2 kg -1 . In red, the soot concentration as a function of albedo required to successfully convert the MODIS albedo field into an impurity field. For both lines, clear-sky conditions are assumed for a SZA of 60 o and RACMO2 irradiance profiles are used to convert narrowband to broadband albedo. (d) Spectral albedo for clean fresh snow (in black), and for an ice profile with the fitted SSA of 0.788 m 2 kg -1 for various impurity concentrations. The first twelve spectral bands of RACMO2 are indicated by vertical dotted lines and black numbers. Red bars and numbers indicate the seven MODIS spectral bands. The albedo for the cases with soot concentrations of 0.2 and 1.5 mu g g -1 are indicated with corresponding colored dots in (c).
Page 5: Hence, albedos ranging from 0.30 to 0.55 observed for the GrIS are obtained by increasing the soot content, with the absorption cross section for soot that is determined by Kokhanovsky (2004) Figure 2. Put labels on the latitude & longitude lines (also on Figure 1). As requested, we added latitude and longitude labels on all Figures with maps.