Impact of updated radiative transfer scheme in RACMO2.3p3 on the surface mass and energy budget of the Greenland ice sheet

This study evaluates the impact of a new snow and ice albedo and radiative transfer scheme on the surface mass and energy budget for the Greenland ice sheet in the latest version of the regional climate model RACMO2, version 2.3p3. We also evaluate the modeled (sub)surface temperature and snow melt, as subsurface heating by radiation penetration now occurs. The results are compared to the previous model version and are evaluated against stake measurements and automatic weather station data of the K-transect and PROMICE projects. In addition, subsurface snow temperature profiles are compared 5 at the K-transect, Summit and southeast Greenland. The surface mass balance is in good agreement with observations, and only changes considerably with respect to the previous RACMO2 version around the ice margins and in the percolation zone. Snow melt and refreezing, on the other hand, are changed more substantially in various regions due to the changed albedo representation, subsurface energy absorption and melt water percolation. Internal heating leads to considerably higher snow temperatures in summer, in agreement with observations, and introduces a shallow layer of subsurface melt. 10


Introduction
The Greenland ice sheet (GrIS) has been losing mass at an accelerating pace in the last decade (Box and Colgan, 2013;Kjeldsen et al., 2015;Bevis et al., 2019;Shepherd et al., 2020). Both surface runoff and solid ice discharge have increased, enhancing mass loss (Bigg et al., 2014). The relative contribution of surface processes with respect to ice discharge, however, has recently increased considerably (Enderlin et al., 2014;Kjeldsen et al., 2015;Mouginot et al., 2019), augmenting the need to accurately 15 model the surface mass balance (SMB) in regional and global climate models.
The SMB, which is the difference between precipitation and ablation, i.e., runoff, sublimation and drifting snow erosion, is highly variable in space and time. Snow and ice melt typically dominate the SMB around the margins of the GrIS, leading to extensive runoff and mass loss of up to 3 m water equivalent (w.e.) yr −1 , while snow fall dominates in the interior (Van den Broeke et al., 2016). In the interior, melt events do not necessarily lead to runoff, as a significant fraction of all melt water 20 refreezes locally (Steger et al., 2017). The length and spatial extent of melt events vary from year to year, but record-breaking melt events like during the summer of 2012 and 2019 can bring melt to almost the entire ice sheet (Nghiem et al., 2012;Bennartz et al., 2013;Tedesco et al., 2013;Tedesco and Fettweis, 2020), inducing a darkening of the snow pack and a change in snow structure. As similar melt events are expected to become more common in a warmer climate (IPCC, 2013;Shepherd et al., 2020), it is therefore important to adequately resolve the individual SMB components. In-situ observations, however, lack tute for Marine and Atmospheric Research Utrecht (IMAU), is adjusted for glaciated areas by introducing a dedicated glaciated tile that includes snow and ice processes and more complete ice-atmosphere interaction (Noël et al., 2015).
Two major components have been updated in the new version Rp3 compared to the previous version Rp2: the multilayer firn module and the snow and ice albedo parameterization. More specifically, the multilayer firn module has increased the vertical resolution and the merging routine reduces mixing of layers with distinct characteristics, reducing numerical diffusion.

70
Rp3 now typically has 50 to 60 active snow layers, with 100 as a maximum. Model output is only available for the first 20 layers and every 3 hours. The plane-parallel broadband snow albedo scheme based on Gardner and Sharp (2010) is replaced by the Two-streAm Radiative TransfEr in Snow Model (TARTES, Libois et al., 2013), which is coupled to RACMO2 with the Spectral-to-NarrOWBand ALbedo (SNOWBAL) module version 1.2 (Van Dalum et al., 2019). TARTES uses the asymptotic radiative transfer theory (Kokhanovsky, 2004;He and Flanner, 2020) and the radiative transfer equation (Jiménez-Aquino 75 and Varela, 2005) to calculate a spectral albedo and subsurface energy absorption by using the geometric-optics method.
SNOWBAL has been developed to couple the output of TARTES with the 14 contiguous shortwave spectral bands of the IFS physics scheme embedded in RACMO2, taking into account sub-band variations of both the albedo and the irradiance.
Bands 13 and 14 are, however, excluded from calculations, as the albedo for these bands can be safely assumed to be zero (Gardner and Sharp, 2010) and all energy contributes to the SEB. Additionally, a new bare ice albedo parameterization using 80 both TARTES and SNOWBAL has been developed. A more detailed model description and evaluation of the snow and ice albedo product can be found in Van Dalum et al. (2020).

Surface mass balance and surface energy budget
The surface energy budget (SEB) of RACMO2, with fluxes toward the surface defined as positive, is defined as 85 with M the surface melt flux (M = 0 when surface temperature T s < 273.15 K), LW d , LW u , SW d and SW u the downward and upward longwave radiative fluxes and downward and upward shortwave radiative fluxes, respectively, SHF and LHF the sensible and latent heat fluxes and G s the subsurface conductive heat flux, all in W m −2 . Melt water is allowed to percolate into deeper layers using the tipping-bucket method, i.e., water fills a layer until the irreducible water saturation is reached (Coléou and Lesaffre, 1998). Excess water moves to the next unsaturated layer, where it can either refreeze, be retained or runs off.
Here, we adopt the following definition of the surface mass balance (SMB), in mm w.e. yr −1 : with SN snowfall, RA rain, SU sublimation, ER drifting snow erosion and RU runoff. RU represents rain and meltwater that 95 is not refrozen or retained in the firn layer. Strictly speaking, this definition of SMB represents the climatic mass balance, as internal accumulation is included (Cogley et al., 2011).

Internal energy absorption
In Rp3, the albedo and energy absorption profiles within snow layers are computed on full-radiation (FR) time steps, which is every whole hour. For all other time steps, the albedo and absorption profiles of the previous FR time step are used as long as 100 the sun is above the horizon. The internal energy absorption that is calculated by TARTES, however, is only valid on a FR time step, as snow layers and radiation entering the snow pack change within the hour.
To be able to calculate internal energy absorption at non-FR time steps, we assume that the net downward shortwave energy flux F (z) decays exponentially within a model layer with attenuation length τ , where depth z > z top , with z top the depth of the upper interface and F top the net downward shortwave flux at the top of the snow layer. F (z) is then: At each FR time step, τ is determined for each snow layer using Eq. (3) and the modeled absorption profile of TARTES. During non-FR time steps, these τ 's are used to distribute the net absorbed shortwave radiation over the model layers. This procedure is repeated for all spectral bands of Rp3, as both the albedo and τ depend on wavelength. For visible light, only a small fraction 110 of incoming radiation is absorbed in the snow pack, but penetrates relatively deeply. For near-infrared (near-IR, 750-1400 nm) radiation, a large fraction is absorbed, but penetration and absorption is limited to the upper layers (Ebert et al., 1995;Gardner and Sharp, 2010;He et al., 2017He et al., , 2018He and Flanner, 2020).
Due to the absorption behaviour for near-IR radiation, a major fraction of incoming radiation is absorbed in the upper millimeters of the snow pack. As for this length scale, the time scale for heat diffusion to the surface is shorter than the model 115 time step (4 minutes) and near surface shortwave radiation absorption therefore needs to be accounted for in the SEB. In order to distinguish between surface and internal energy absorption, we assume that the fraction of absorbed energy attributed to the SEB is 1 at the surface and linearly decreases to 0 at the maximum skin layer equilibration depth (SLED). Using scale analysis, we estimate the SLED for both snow and ice to be 5 mm (Appendix A). All energy absorbed beyond the SLED can The absorbed energy that contributes to the SEB is the difference between the total energy absorbed in this layer and E internal .

125
For the layer in which the SLED is located, Eq. (4) is evaluated to z sled ; all energy below SLED is counted as internal energy absorption. Figure 1 illustrates F (z), dF (z)/dz and the energy contributing to E SEB and E internal as a function of depth for a typical fresh snow (Fig. 1a) and ice layer (Fig. 1b) for Rp3's spectral band 6 (778-1242 nm) and 4 (442-625 nm), respectively. As energy absorption depends on F (z) and τ , which in turn depends strongly on wavelength and snow structure (Meirold-Mautner 130 and Lehning, 2004; Ackermann et al., 2006;Warren et al., 2006;He et al., 2017He et al., , 2018He and Flanner, 2020;Cooper et al., 2020), significantly more energy is absorbed within the SLED for fresh snow and near-IR radiation, as τ is small for this case, than for ice and visible light, for which τ is large, resulting in more internal energy absorption close to surface and a larger E SEB .

135
In this study, we test Rp3 on a 11 km grid for the GrIS and its immediate surroundings between 2006 and 2015, using September 2000 to December 2005 as spinup. At the lateral boundaries, Rp3 is forced with ERA-Interim data (Dee et al., 2011). The only impurity type considered is soot, which is homogeneously distributed for all snow and firn layers. The impact of soot is known to be underestimated, so a fixed concentration of 5 ng g −1 is used, which is higher than the typical 3 ng g −1 soot concentration of the interior (Chylek et al., 1992;Doherty et al., 2010;Dang et al., 2015). A spatially variable soot concentration is used for The firn-column, i.e, snow and ice density, thickness, temperature, grain radius and water concentration, is initialized for all active layers with output of Rp2 (Noël et al., 2018b), and all grains are spherically shaped. Bare ice is identified if the continuous set of layers, counted from bottom up, have a density larger or equal to 899 kg m −3 .

145
In order to highlight the impact of model changes, we investigate the sensitivity of various parameters by comparing the results with a run of Rp3 that is without internal energy absorption (Rp3 WIE). Rp3 WIE also covers the same period of analysis as Rp3 and Rp2. Furthermore, we investigate the sensitivity of the numerical choices in the implementation of internal heating.   between -500 and 500 mm w.e. yr −1 and 500 mm w.e. yr −1 elsewhere. (b) SMB difference (Rp3 -Rp2). A linear color scale with a step size of 100 mm w.e. yr −1 between -500 and -100 and between 100 and 500 mm w.e. yr −1 is used and 20 mm w.e. yr −1 elsewhere. The dots in the southwest are the locations of the K-transect AWS stations, from left to right: S6, S9 and S10. KAN-U and KAN-M of the PROMICE data set are located close to S10 and between S6 and S9, respectively, but are not shown separately. Subsurface temperature profiles are available for Summit. The star indicates a location that will be discussed in Sect. 6.1. (c) The relative difference between Rp3 and Rp2. Patterns A and B are discussed in the text.

150
Three type of measurements are used in this study: stake measurements to determine the SMB, AWS observations to evaluate the SEB, 2-m temperature (T 2m ) and 10-m wind speed (v 10m ), and subsurface temperature observations. The SMB observa- For the evaluation of SMB and SEB, the two closest grid points to an observational site are selected in RACMO2 and linearly 170 interpolated between them. Interpolation between two grid points is, however, not desirable for the analysis of subsurface processes, as the size and depth of snow layers in RACMO2 can alter considerably between two neighboring points. Instead, we use the nearest model grid point for subsurface processes.

180
Integrated over the ice sheet, the SMB increases by 17.0 mm w.e. yr −1 , or 9.8%. In the interior, significant differences are observed only in south Greenland, with a considerable increase of snow melt (Fig. 3a). As radiation penetration now allows internal heating, it consequently raises subsurface temperatures and increases melt. Because all snow melt refreezes in the snow pack, these changes do not affect the local SMB ( Fig. 3b and c). Refreezing, however, changes the structure of the snow pack, inducing more energy absorption and higher snow temperatures, leading to more melt.

185
The outer rim of the ice sheet, except in the southeast, is characterized by a strong SMB increase. In Rp2 at 11 km, the bare ice albedo of the majority of the outermost glaciated points is 0.30 due to contamination with tundra albedo, causing too much melt and runoff. This artifact is mitigated for higher resolutions and is solved in Rp3 (Van Dalum et al., 2020), lowering the snow melt and runoff and increasing the SMB. In the southeast close to the ice-sheet margin, this issue is absent, as at this low resolution the number of grid points for which bare ice reaches the surface is limited. Snow melt is slightly higher and 190 refreezing is similar to Rp2 for this region, resulting in more runoff. This is discussed in more detail in Sect. 5.  The lower ablation area of the GrIS in the southwest experiences a strong increase in refreezing (Fig. 3b) some regions in the southwest show a small increase in snow melt, which is induced by, on average, a slightly lower albedo in Rp3 (Van Dalum et al., 2020). As a result, the increase in snow melt and refreezing balance out, leading to a runoff similar to Rp2 (Fig. 3) and little change in SMB (Fig. 2b).
In Rp3, snow melt has increased in the percolation zone and most of the additional melt water refreezes locally. Close to the equilibrium line, however, the melt water buffering capacity is exceeded due to the addition of melt water and runoff occurs, lowering the SMB (pattern A in Fig. 2c). Similarly, in the upper ablation zone, where the winter snow cover lasts long during summer, enhanced melt leads to a small reduction of refreezing, as the refreezing capacity is consumed faster.
The northeast shows a large area with reduced SMB (pattern B in Fig. 2c). This region is characterized by dry conditions, resulting in a thin winter snow layer that melts away quickly in spring, exposing bare ice. Due to the low albedo of bare ice and radiation penetration, the ice temperature is raised, leading to melt in multiple model layers reaching greater depths and 205 prolonging the melt season (discussed in more detail in Sect. 6.1), consequently enhancing melt and runoff (Fig. 3a, c, d and f). Furthermore, it usually takes several weeks after the ablation season ends for fresh snow layers to grow to a thickness of over 10 cm. During this time, some radiation can still reach bare ice layers, lowering the albedo and inducing more energy absorption.
To summarize, averaged over the ice sheet, snow melt has barely changed, i.e., -1.1 mm w.e. yr −1 (-0.33%), compensating 210 the increase of melt in south Greenland with the reduced melt at the margins. The runoff, on the other hand, is lowered by -18.7 mm w.e. yr −1 (-8.6%) while refreezing has increased by 18.1 mm w.e. yr −1 (12.7%). By percentage, snow melt and refreezing have increased considerably in the interior of the ice sheet, except for a large area in the center (Fig. 3d and e). This is discussed in more detail in Sect. 3.2. No runoff is modeled in the interior ( Fig. 3c and f), as melt water refreezes locally.
Around the K-transect, the relatively large increase of refreezing consequently lowers the modeled runoff. A relatively strong and root-mean-square error (RMSE, in m w.e. yr −1 ) are also displayed. Results for S6 are shown with squares.

Cloud cover and melt
For a large area of central GrIS, snow melt and refreezing have decreased considerably by percentage ( Fig. 3d and e). Melt events in this region, although rare, are almost always associated with above-average cloud cover (Fig. 4a), reducing the amount of SW radiation that reaches the surface. In general, the albedo of Rp3 is lower than Rp2 for this region (Van Dalum et al., 220 2020), but for cloudy conditions and during melt events, the albedo of Rp3 is considerably higher (Fig. 4b), leading to less energy absorbed in the snow pack and subsequently less melt and refreezing. This albedo increase is mostly caused by the changing spectral distribution of irradiance during cloudy conditions, as relatively more irradiance is ultraviolet (UV) and visible light, for which the albedo is high, and relatively less is IR radiation, for which the albedo is low (Van Dalum et al., 2019). A large fraction of UV and visible radiation is absorbed internally, heating subsurface layers that can potentially increase 225 melt. This is discussed in more detail in Sect. 5.2. This effect, however, is not enough to mitigate the reduced irradiance and higher albedo and therefore less snow melt is modeled. Note that the total amount of snow melt modeled for this region in central GrIS is on average very small (2 mm w.e. yr −1 ) compared to ablation areas (e.g., 2500 mm w.e. yr −1 at S6). The black line is the 1-on-1 line and the red line shows the orthogonal total least squares regression of the data, with b0 the slope and b1 the intercept. Number of records (N), correlation coefficient (R 2 ), root-mean-square error (RMSE, in W m −2 ) and bias are also displayed. In all panels, positive values imply an energy flux directed towards the surface. Figure 5 shows the SMB of Rp3 and Rp2 with respect to various SMB observations. The correlation coefficient, root-mean-230 square error (RMSE) and bias with respect to observations are also shown. A slightly higher SMB is modeled for Rp3 compared to Rp2, reducing the bias with observations. The differences with observations are particularly relevant for the measurement sites in the lower ablation zone (S6, squares in Fig. 5) and are caused by increased refreezing in Rp3 (Fig. 3b). Note that the spread is considerably higher for the ablation zone, illustrating that the large temporal variability is not always captured properly with this resolution in RACMO2. Increasing the resolution might improve the modeled estimates for the locations 235 close to the ice margin.   (Table 1). The small bias of LWd 240 and shortwave downward radiative flux (SWd) shows that the influence of cloud cover is captured well. There is a tendency for Rp3 to underestimate LWd during cold and dry, cloud-free winter (small SWd) conditions, resulting in underestimated LWu for the same cold days.

SMB Observations
SWd and upward shortwave radiative flux (SWu) are both in good agreement with observations. The R 2 , bias and RMSE of these fluxes have improved in Rp3 compared to Rp2 (Table 1), which confirms that the albedo has improved, in agreement Both the sensible and latent heat flux (SHF and LHF, respectively), show a considerable spread ( Fig. 6c and f). While the LHF is small and contributes little to the SEB, the SHF is important to model melt events properly. The difference with observations for the 2-m temperature (T 2m ) 2-m and 10-m wind speed (V 10m ) ( Table 1 and Fig. B1), however, are relatively small, and can therefore not explain the large spread of the turbulent fluxes. As is discussed in Noël et al. (2018b), the bias in 250 SHF can be attributed to an inaccurate representation of the roughness length, especially for bare-ice conditions. Figure 7 shows the climatology of the SEB components of in-situ, Rp2 and Rp3 data for S6, S9 and S10. The Rp3 and Rp2 radiative and turbulent fluxes are generally similar, and agree with in-situ observations. Some differences with observations are, however, still observed. The net absorbed longwave radiation (LWn), defined as LWd -LWu, is significantly lower during winter for both Rp2 and Rp3, which is mostly due to the temperature error induced by an overestimation of SHF and partly  results as Rp2 are therefore expected (Noël et al., 2018b). Improving the representation of surface roughness may also improve the results. During summer, only the net absorbed shortwave radiation (SWn), defined as SWd -SWu, shows considerable discrepancies for S6 and S9. The introduction of a new radiation penetration and snow albedo scheme reduces the differences with observations for Rp3 compared to Rp2. This is especially noticeable at the onset of the accumulation season at S6 when 260 a thin fresh snow layer forms on top of bare ice. As Rp2 does not include radiation penetration, the ice layers are colder and the melt season ends earlier, allowing the early formation of a fresh snow layer, raising the albedo. The albedo differences are relatively small for S10, as this location is characterized by a thick firn layer. The evaluation of the snow albedo scheme of Rp3 is discussed in more detail in Van Dalum et al. (2020).

K-transect
Internal energy absorption of shortwave radiation is an important addition in Rp3. We illustrate the absorption of solar radiation and grain radius as a function of depth for the K-transect for 2012 in Fig. 8. S6 is situated in the ablation zone and this dry region is characterized by a relatively thin fresh snow pack after winter. During June, the ice horizon quickly moves up exposing bare ice until September. As S6 is situated in the so-called 'dark zone' (Van de Wal and Oerlemans, 1994;Wientjes et al., 2011), 270 the exposed bare ice has a high impurity concentration and a large grain radius, leading to extensive internal energy absorption up to 15 cm deep. During the accumulation season, internal energy absorption is limited, except during sporadic melt events when the surface grain radius grows rapidly, allowing radiation to penetrate to deeper layers.
At S10, situated in the accumulation zone, higher winter accumulation provides a thicker snow pack at the start of the melt season than at S6. Consequently, the whole snow column only melts away during extreme melt events like in the summer 275 of 2012 and 2019, in which bare ice surfaces for a few weeks in August (Fig. 8b). In late June and July, firn layers that are characterized by a large grain radius reach the surface and induce significant absorption of solar radiation, although less than for bare-ice conditions. The absorption, however, is limited to the upper 5 cm. Note that the formation of a thin fresh snow layer in early July considerably diminishes internal energy absorption.

Distribution of energy 280
Since a significant fraction of the shortwave radiation is absorbed in the upper few millimeters of the snow column, not all incoming radiation contributes to internal heating (Sect. 2.3). The fraction of shortwave radiation absorbed that directly contributes to the SEB depends on wavelength and snow conditions, which is illustrated in Fig. 9. This figure shows, as an example, the absorption of energy for a summer day for a grid point at Summit in central Greenland (left bars) and S6 (hatched right bars) for the first twelve spectral bands. Band 13 and 14 are not shown, as all energy for these bands is added to the SEB.

285
Blue shows the energy absorbed internally and orange the energy that contributes to the SEB, with numbers above the bars indicating the percentage of surface radiation that is absorbed. At Summit, almost all radiation is absorbed internally for visible light (bands 1-5), but the amount of radiation that remains in the snow pack for these bands is limited. For IR radiation (bands 6-12), considerably more absorption takes place. The attenuation length in snow is short for IR radiation (He and Flanner, 2020), leading to a strong absorption within the SLED, 290 increasing the fraction of energy that contributes to the SEB.
S6 has bare ice at the surface for the selected day, and considerably more radiation is absorbed in the visible light bands due to a high concentration of soot. Most notably, is the large increase of energy absorbed in band 6. This is a wide band with a considerable amount of energy available. It is especially sensitive to albedo drops induced by an increase of grain radius and density, e.g., when bare ice reaches the surface (Van Dalum et al., 2019).

295
For bands 7 to 12, almost all incoming radiation is absorbed at S6 due to the large grain radius and density. For these bands a larger fraction of energy contributes to the SEB at S6 compared to Summit. As the grain size at S6 has increased more strongly than the density when compared to the snow pack at Summit, the number of grains in the SLED is therefore lower. In other words, the snow pack at S6 misses a scattering fresh snow top layer, hence increasing the chance of near-IR radiation to penetrate through the SLED without being absorbed. This effect is less pronounced for visible light, as the SLED is so small 300 that visible light travels in and out virtually without loss, especially if the soot concentration is as low as that at Summit.

Sensitivity experiments
In this section, we discuss the impact of internal energy absorption on the subsurface temperature for various locations in winter and summer. Furthermore, we investigate the 10-m snow temperature (T 10m ) and the impact of internal energy absorption on the SMB by comparing the results with Rp3 without internal energy absorption (WIE).

Subsurface temperature
First, we analyze the impact of internal energy absorption on subsurface snow temperatures. Figure 10 shows the subsurface temperature profile for Rp2, Rp3 and Rp3 WIE for S6 and Summit for a typical winter day and summer day, and for an extraordinary warm summer day for S6 and a grid point in southeast Greenland (locations are shown in Fig. 2b). As the modeled temperature is not interpolated between layer mid-points, the stepwise temperature profile indicates individual model 310 layers. Rp3 generally has thinner snow layers, leading to a higher vertical resolution, which is especially relevant close to the surface. As Fig. 10 only shows the upper 20 snow layers, the temperature curve of Rp3 therefore does not reach as deep as Rp2.
17 https://doi.org/10.5194/tc-2020-259 Preprint. Discussion started: 30 September 2020 c Author(s) 2020. CC BY 4.0 License. During winter (Fig. 10a), the temperature profiles are similar for both S6 and Summit. Internal heating is limited, as only a small amount of radiation reaches the surface and both sites are covered with fine-grained fresh snow layers. Only a significant 315 change is observed in the upper first centimeters, with thinner snow layers and higher temperatures in Rp3.
During summer (Fig. 10b), internal energy absorption has a clear impact on subsurface temperatures. For S6, melting occurs in multiple layers reaching to a depth of about 14 cm, which coincides approximately with depths that absorption of solar radiation is still relevant (Fig. 8c). With the addition of internal heating, the temperature of deeper layers is raised as internal heat cannot dissipate as easily as on the surface. Higher subsurface snow temperatures make the snowpack more susceptible to 320 internal melt, thus increasing the vertical melt extent. Note that all melting layers are still connected to the surface, effectively reducing it to a single melt layer. The impact of internal heating is also illustrated by the similar temperature profile of Rp3 WIE with respect to Rp2.
For Summit, much colder conditions prevail during summer. Compared to Rp2, the subsurface temperatures of Rp3 are considerably higher (up to 5 • C) and match the observations better, but are slightly overestimated in the upper 10 centimeters. Rp3 WIE shows lower temperatures similar to Rp2, leading to a cold bias with the observations. The temperature decline with depth, however, is more gradual for Rp3 WIE compared to Rp2. Adding internal energy absorption considerably increases the ability of RACMO2 to reproduce realistic subsurface temperature profiles for Summit. Figure 10c shows the subsurface temperature profile during a melt event in the accumulation zone of southeast Greenland (SE GRL, star in Fig. 2b)  3d and e, which shows a relatively strong increase in snow melt and refreezing around central GrIS for Rp3 compared to Rp2, especially in the southeast where this grid point is located. As this location in SE GRL is too far from the ice margin, no runoff is modeled and all melt water refreezes locally. Even though the snow melt and consequent refreezing are still small and do 340 not alter the SMB, they do change the snow structure.
We also investigate the sensitivity of the subsurface snow temperature to the the skin layer equilibration depth z sled , which we introduced in Sect. 2.3 and set to 5 mm using scale analysis. Underestimating z sled leads to an overestimation of internal heating, as not enough heat diffuses to the surface within the model time step, trapping too much heat beneath the surface.
Increasing z sled only has a limited effect on the subsurface temperature (Appendix Fig. B2). As expected, lowering z sled 345 results in too much internal heating and raises the snow temperature, illustrating that an underestimation of z sled should be avoided.

Snow temperature at 10 m depth
In the absence of melt and refreezing, the T 10m is regarded as a good proxy of the mean annual surface temperature. Here, we analyze whether this assumption is true given the observed impact of radiation penetration on summer subsurface temperatures.
350 Figure 11 shows the averaged yearly mean T 10m at the end of the simulation for Rp3 (Fig. 11a), the difference between the average skin temperature (T skin ) and T 10m (Fig. 11b) and the T skin − T 10m difference for Rp3 and Rp3 WIE with respect to Rp2 ( Fig. 11c and d, respectively).
In the interior, the difference between T skin and T 10m is small for Rp3, illustrating that T 10m is indeed a good proxy of the mean annual surface temperature for this region. For Rp3, T skin −T 10m is comparable to Rp2 (Fig. 11c), for Rp3 WIE, however,

355
it is considerably larger (Fig. 11d). As all solar energy is absorbed at the surface for Rp3 WIE and the albedo is slightly lower compared to Rp2, T skin is consequently increased while the energy that reaches 10 m depth is decreased, thus lowering T 10m .
In other words, Fig. 11d highlights the importance of radiation penetration. In the percolation zone (A of Fig. 2c), the difference between T skin and T 10m is considerable (Fig. 11b), owing to latent heat being released upon refreezing. Consequently, T 10m is not representative for the surface temperature in this zone. There are 360 also large differences in southeast Greenland, which is characterized by firn aquifers that keep the deeper snow temperature near the melting point year round.
To conclude, T 10m is a good measure for the surface temperature, except close to the ice-sheet margin in the southeast and the percolation zone. Moreover, subsurface energy absorption reduces T skin − T 10m in the interior.

365
In this subsection, we discuss the impact of radiation penetration on the SMB and analyze its components by comparing Rp3 WIE with Rp3. Figure 12 shows the SMB, snow melt and refreezing difference between Rp3 WIE and Rp3. Note that in Rp3 WIE radiation penetration is still used to determine the albedo, but all absorbed shortwave radiation is added to the SEB.
The SMB difference of Fig. 12a is limited in the interior, where melt does not run off. The SMB is lower for Rp3 WIE in most of the ablation zone (Fig. 12a), in particular in the southwest (area C of Fig. 12c), where we observe a relatively small 370 snow melt increase (Fig. 12b). More importantly, however, is that internal radiative heating reduces refreezing in the ablation zone (Fig. 12c), especially in the southwest, but also in the northwest and northeast (area C, D and B of Fig. 2c, respectively), so that runoff increases, reducing the SMB. The SMB difference of Fig. 12a is limited in the interior, where melt does not run off.
In the central part of the interior, the melt and refreezing differences are limited, supporting the findings of Sect. 3.2. In a broad elevation band around the center, less snow melt and refreezing is observed, illustrating the importance of subsurface energy 375 absorption, which prevents the surface from reaching the melting point. In south Greenland, the melt difference is considerably smaller than the difference between Rp3 and Rp2 (Fig. 3d), although still not negligible, illustrating that internal heating is not the only relevant process here. This region is also characterized by an albedo decrease (Van Dalum et al., 2020), resulting in more absorption of shortwave radiation and snow melt, contributing significantly to the observed difference between Rp3 and Rp2.

Conclusions
In this study, we evaluated the impact of a new snow and ice albedo and radiative transfer scheme that allows internal heating on the surface mass balance (SMB), surface energy budget (SEB) and firn and snow temperature of the GrIS in the latest version of the regional climate model RACMO2. In this new version Rp3, we also updated the multilayer firn module (Van Dalum et al., 2020) and developed a method to distribute total absorbed shortwave radiation between the SEB and to internal absorption. 385 We have run Rp3 for the GrIS for 2006-2015 using ERA-Interim data at the lateral boundaries. The modeled SMB correlates well with ablation stake measurements and the SEB is in good agreement with AWS observations at the K-transect and for a selection of PROMICE stations. Moreover, subsurface radiative energy absorption leads to higher summer snow temperatures and a larger vertical melt extent, improving the agreement with observed snow temperatures at Summit. Compared to the previous RACMO2 version Rp2, the SMB difference is small in the interior, as any change in snow melt is 390 balanced by refreezing. In the percolation zone close to the equilibrium line, the SMB is lowered as the melt water buffering capacity is exceeded due to the addition of melt, leading to runoff. In the ablation zone in the southwest, considerably more refreezing occurs due to pore space created in ice by melt induced by subsurface energy absorption. Finally, the SMB is higher at the very margins due to the removal of an artifact in the ice albedo of Rp2, which reduces snow melt. These differences are small or absent in the wet southeast, where bare ice rarely reaches the surface. For this region, relative differences in snow melt 395 and refreezing are substantial due to subsurface energy absorption and melt water percolation. This increased snow melt and refreezing does not lead to significant SMB changes, but does change the snow structure. Furthermore, subsurface temperature profiles show a more gradual temperature decrease with depth in Rp3 during summer in the interior, and melt reaching deeper layers both in the accumulation and ablation zone.
There is, however, still room for improvement. The resolution used in this study, i.e., 11 km horizontal resolution, is sufficient 400 to resolve the SMB and SEB in the interior of the ice sheet, but is insufficient to resolve the inhomogeneous ablation zone, especially close to the ice margin. For example, bare rock cropping out of the ice, steep slopes or snow patches are usually smaller than the current RACMO2 resolution. More fundamentally, the ice albedo, although significantly improved in Rp3, could be further improved upon by using a dynamic Mie-scattering theory model (Gardner and Sharp, 2010). In addition, a new impurity scheme that would allow the introduction of more impurity types and varying concentrations for each layer in 405 space and time and would be beneficial.
To conclude, Rp3 is capable to accurately reproduce the SMB, SEB and (sub)surface temperature when compared to insitu observations. Differences with Rp2 are most pronounced in the ablation zone close to the ice margin, in south Greenland and in the percolation zone. Furthermore, the snow temperature with depth is altered ice-sheet wide due to subsurface energy