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

Radiative transfer in snow and ice is often not modeled explicitly in regional climate models. In this study, we evaluate a new englacial radiative transfer scheme and assess 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 melt, as radiation penetration now enables internal heating. 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 5 addition, subsurface snow temperature profiles are compared at the K-transect, Summit and southeast Greenland. The surface mass balance is in good agreement with observations, with a mean bias of -31 mm w.e. yr−1 (-2.67%), and only changes considerably with respect to the previous RACMO2 version around the ice margins and near the percolation zone. 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 higher snow temperatures in summer, in 10 agreement with observations, and introduces a shallow layer of subsurface melt. Hence, this study shows the consequences and necessity of radiative transfer in snow and ice for regional climate modeling of the Greenland ice sheet.


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 15 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 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 leading to extensive runoff typically dominates the SMB around the 20 margins of the GrIS, leading to mass loss of up to 3 m water equivalent (w.e.) yr −1 , while snow fall dominates in the interior . In the interior, melt events do not necessarily lead to runoff, as a significant fraction of all melt water refreezes locally (Steger et al., 2017). The length and spatial extent of melt events vary from year to year, but recordbreaking melt events like during the summer of 2012 and 2019 can bring melt to almost the entire ice sheet (Nghiem et al., SEB for the GrIS in Rp3. Section 2 discusses the model and initialization, expands on the concepts of SMB, SEB and internal energy absorption, and discusses the in-situ data sets in more detail. In Section 3, internal energy absorption due to radiation penetration is assessed. In Section 4, sensitivity experiments further highlight the impact of internal heating on the subsurface temperature. Results are compared with observations and compared to the previous RACMO2 version, 2.3p2 (Rp2) (Noël et al., 2018b) and to a Rp3 model version without internal energy absorption (Rp3 WIE). In Section 5, results are analyzed by evaluating the SEB to observations and by comparing modeled SEB of Rp3 with Rp2. Similar to Section 5, Section 6 focuses on the SMB and its components. We also compare the SMB of Rp3 with Rp3 WIE. Finally, Section 7 summarizes the results and conclusions are drawn.
2 Methods and data 2.1 Regional climate model  Undén et al., 2002). The polar (p) version of RACMO2, which is developed at the Institute for Marine and Atmospheric Research Utrecht (IMAU), is adjusted for glaciated areas by introducing a dedicated glaciated 75 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, which makes use of a radiative transfer scheme.

Multilayer firn module
Four modifications have been applied to the multilayer firn module, which we briefly address here and are discussed in more 80 detail in Van Dalum et al. (2020). Firstly, Rp2 uses a prognostic fresh snow layer, which is effectively a sublayer of the uppermost snow layers. This sublayer is removed, and the uppermost snow layers are now allowed to be as thin as a few millimeters. Secondly, if merging is necessary, a snow layer now merges with its most similar adjacent layer instead of the next layer. Numerical diffusion is prevented by redistributing mass if a thin layer merges with a thick layer. Layers formed by local accumulation are not allowed to merge with glacial ice. Furthermore, the effective vertical resolution has increased, with 85 typically 50 to 60 active snow layers, up to a maximum of 100. Model output, however, is only available for the first 20 layers.
Thirdly, internal melt will thin a subsurface snow layer if the density is lower than 700 kg m −3 . Pore space is created for ice layers with a density of more than 830 kg m −3 . A combination of both occurs for firn with intermediate densities. Melting of the upper layer always leads to thinning. Finally, the initialized ice density is changed from 910 in Rp2 to 917 kg m −3 in Rp3.
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). The broadband albedo parameterization of Gardner and Sharp (2010) is based on tuning parameters and lookup tables; and parameterizes the albedo impact of SZA, wavelength of irradiance, grain radius, cloud cover and impurities. Only the upper snow layer is considered for its calculations, with the 95 second layer as a semi-infinite background layer. As the upper snow layers in RACMO2 are often only millimeters thick, there is virtually no radiation penetration in Rp2.
The albedo scheme of Rp3 is fundamentally different, as TARTES uses the asymptotic radiative transfer theory (Kokhanovsky, 2004;He and Flanner, 2020) and the radiative transfer equation (Jiménez-Aquino and Varela, 2005) to calculate a spectral albedo and subsurface energy absorption by using the geometric-optics method. Spectral radiative transfer in TARTES de-100 pends on snow layer density, impurity concentration, grain radius and grain shape. Here, grains are assumed to be spherically shaped; and prognostic estimates are provided for the other variables by the multilayer firn model. 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. SNOWBAL selects a predefined representative wavelength for given atmospheric conditions suitable for TARTES to use, such that a narrowband albedo and 105 subsurface energy flux can be accurately represented by one spectral evaluation of TARTES (Van Dalum et al., 2019). The representative wavelength depends on the SZA and vertically integrated water vapor for clear-sky conditions and ice and liquid water path for cloudy conditions. Bands 13 and 14, covering 3077-3845 and 3846-12500 nm, respectively, are 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.

Ice albedo
Additionally, a new bare ice albedo parameterization using both TARTES and SNOWBAL has been developed. In Rp2, the bare ice albedo is prescribed by the lowest 5% of the 16 day diffuse albedo product of 1 km MODIS data (MCD43A3v5, Schaaf and Wang, 2015), resampled to the RACMO2 grid and limited between 0.3 for dark ice and 0.55 for perennial snow (Noël et al., 2018b). In Rp3, we replace these predefined bare ice albedos by using TARTES and SNOWBAL for each spectral 115 band, allowing for varying bare ice albedo and estimates of subsurface heating. TARTES, however, is not suitable to be applied directly to ice, as it is not based on Mie-scattering theory, and some approximations have to be made.
Firstly, we determine the specific surface area (SSA) of a semi-infinite ice layer such that TARTES and SNOWBAL calculates a broadband albedo of 0.6, which represents the broadband albedo of clean blue ice (Reijmer et al., 2001;Dadic et al., 2013) for clear-sky conditions and a reference SZA of 60 • . For this, we find an SSA of 0.788 m 2 kg −1 (4.152 mm grain size).

120
Then, using this SSA and SZA, a broadband albedo can be determined for a range of impurities. With this impurity range, the MCD43A3v5 MODIS albedo can be converted into an impurity concentration for each model grid point. Each time a layer is identified as bare ice, TARTES and SNOWBAL use the prescribed SSA and impurity field to do its calculations. For subsurface glacial ice layers, we use this procedure as well. Superimposed ice is treated differently than glacial ice, as it is formed by refreezing of meltwater in snow layers and has a granular structure (Granskog et al., 2006). In Rp3, superimposed ice is 125 therefore treated as a snow layer, with a minimum grain radius of 0.720 mm for a layer density of 750 kg m −3 , and the grain radius increases linearly to the bare ice value of 4.152 mm for a density of 917 kg m −3 . A more detailed model description, evaluation and discussion 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 into deeper layers using the tipping-bucket method, i.e., water fills a layer until 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.
There is no lateral flow between grid points. The subsurface conductive heat flux G s is the energy flux between the skin layer, i.e., the contact layer between the atmosphere and the surface, and the uppermost layer of the firn column. Between all model layers, G s is also derived and used for the subsurface temperature evolution, but this flux is not stored.

140
Here, we adopt the following definition of the surface mass balance (SMB), in mm w.e. yr −1 : with SN snowfall, RA rain, ER drifting snow erosion, SU surface and suspended snow sublimation, and RU runoff. RU represents rain and meltwater that is not refrozen or retained in the firn layer. Strictly speaking, this definition of SMB represents 145 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 of 4 minutes, the albedo and absorption profiles of the previous FR time step are used as long as the sun is above the horizon. The internal energy absorption that is calculated by TARTES, however, is only 150 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 defined positive if downwards and 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 160 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). On average, considerable absorption of solar radiation (more than 20 W m −2 ) is limited to the upper 20 centimeters for clean ice, and only to several milliliters for snow.
Due to the absorption behaviour for near-IR radiation, a major fraction of incoming radiation is absorbed in the upper 165 millimeters of the snow pack. As for this length scale, the time scale for heat diffusion to the surface is shorter than the model time step 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). In other words, the SLED is defined as the maximum depth that some energy can still equilibrate with the surface within a model time step. Using scale analysis,

170
we estimate the SLED for both snow and ice to be 5 mm for the 4 minute time step used in this study (Appendix A). All energy absorbed beyond the SLED can be ascribed to internal absorption. Hence, for a layer ranging from top z top to bottom z bot that is completely above the SLED, i.e., z bottom ≤ z sled , and by using Eq. (3), the absorbed internal energy E internal is given by the dissipated flux: The absorbed energy that contributes to the SEB is the difference between the total energy absorbed in this layer and E internal .
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 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), more energy is absorbed within the SLED for fresh snow and near-IR radiation, as τ is small for this case, than for ice 185 and visible light, for which τ is large, resulting in more internal energy absorption close to surface and a larger E SEB .

RACMO2 simulations
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 190 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 bare ice such that the Rp3 bare ice albedos match the MODIS observations for clear-sky conditions and typical solar zenith angle (Van Dalum et al., 2020).
The firn-column, i.e, snow and ice density, thickness, temperature, optical grain radius and water concentration, is initialized 195 for all active layers with output of Rp2 (Noël et al., 2018b). 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 .
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 and all energy is added to the SEB. Furthermore, we investigate the sensitivity of the numerical 200 choices in the implementation of internal heating by two sensitivity experiments, in which z sled is set to 2.5 and 10 mm.
To determine the statistical significance of the bias between model versions or the bias with observations, we use statistical bootstrapping with a significance of two standard deviations.

In-situ observations
Three type of measurements are used in this study: stake measurements to determine the SMB, AWS observations to evaluate  sensor tilt corrections. For the K-transect, this is discussed in more detail by Smeets et al. (2018). Smeets et al. (2018) also report an uncertainty for shortwave radiation of approximately 1%, for longwave radiation of ±5 W m −2 , for T 2m of max ±0.5 The majority of observations are performed outside the time window of the model run or within two grid points of the 220 ice-sheet margin; these are dismissed, as the ice-sheet margins are not properly captured at the 11 km resolution used for these runs. For example, the modeled surface elevation above sea level for the PROMICE stations NUK-U and THU-U, both located close to the ice margin within two grid points on the RACMO2 grid, show an undesirably large elevation difference with observations: 1069 and 665 m in Rp3 versus 1120 and 760 m observed, respectively. Here, we accept only grid points with an elevation differences of 50 m or less. This limits the data set to S6, S7, S8, S9 and S10 for the K-transect, to  and KAN-U for PROMICE and to two SMB measurement locations slightly to the northwest of S9.
For the evaluation of SMB and SEB, the two closest grid points to an observational site are selected in RACMO2 and linearly 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 change between two neighboring points. Instead, we use the nearest model grid point for subsurface processes. internal energy absorption is limited, except during sporadic melt events when the surface optical grain radius grows rapidly, allowing radiation to penetrate to deeper layers.

240
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 of 2012 and 2019, in which bare ice surfaces for a few weeks in August (Fig. 2b). In late June and July, firn layers that are characterized by a large optical grain radius reach the surface and induce 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 245 in early July diminishes internal energy absorption.  Since a 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. 3. This figure shows, as an example, the 250 absorption of energy for a summer day for a grid point at Summit in central Greenland (left bars) and S6 (hatched right bars)

Distribution of energy
for the first twelve spectral bands. Band 13 and 14 are not shown, as all energy for these bands is added to the SEB. 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 UV and visible light (bands 1-5), but the amount of radiation that 255 remains in the snow pack for these bands is limited. For IR radiation (bands 6-12), 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, increasing the fraction of energy that contributes to the SEB.
S6 has bare ice at the surface for the selected day, and 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 For bands 7 to 12, almost all incoming radiation is absorbed at S6 due to the large optical grain radius and density. Compared to Summit, a smaller fraction of energy contributes to the SEB for these bands. As the grain size has increased more strongly than the density at S6 compared to the snow pack at Summit, the number of grains in the SLED is therefore lower. In other 265 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 that visible light travels in and out virtually without loss, especially if the soot concentration is as low as that at Summit.

Temperature experiments
Compared to observations, the 2-m air temperature difference is only small (Table 1 and Fig. B1). The impact of internal energy 270 absorption on the subsurface temperature, however, is larger for various locations in winter and summer, and is discussed in this section. We also compare the results with Rp3 WIE. Furthermore, we discuss the 10-m snow temperature (T 10m ). We also investigate the impact of z sled on the subsurface temperature, which is shown in Appendix B, Fig. B2.

Subsurface temperature
First, we analyze the impact of internal energy absorption on subsurface snow temperatures. Figure 4 shows the subsurface 275 temperature profiles for Rp2, Rp3 and Rp3 WIE for S6 and Summit for a winter (Fig. 4a) and summer day (Fig. 4b). The figure also shows an example of a melt event in the accumulation zone of southeast Greenland and compares it with a melt event at S6 (Fig. 4c). All locations are shown in Fig. 5c. As the modeled temperature is not interpolated between layer midpoints, the stepwise temperature profiles indicate individual model layers. Rp3 generally has thinner snow layers, leading to a higher vertical resolution, which is especially relevant close to the surface. As Fig. 4 only shows the upper 20 snow layers, the 280 temperature curve of Rp3 therefore does not reach as deep as Rp2. During a typical winter day (Fig. 4a, with a T 2m of -11 and -18 • C respectively in Rp3), 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. We observe only changes in the upper first centimeters, with thinner snow layers and higher temperatures in Rp3.

285
During a typical summer day (Fig. 4b, with a T 2m of 5 and -8 • C, respectively in Rp3), 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. 2c). 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 internal melt, thus increasing the vertical melt extent. Note that all 290 Table 1. Statistics of the daily-averaged Rp2 and Rp3 SEB components, atmospheric 2-m temperature (T2m) and 10-m wind speed (V10m) compared to in-situ observations between 2006 and 2015. Data include S6, S9 and S10 (S10 is not available for T2m and V10m) of the For Summit, 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, 295 however, is more gradual for Rp3 WIE compared to Rp2. Adding internal energy absorption increases the ability of RACMO2 to reproduce realistic subsurface temperature profiles for Summit. snow and firn where melt water can easily percolate, similar to the snow pack illustrated for S10 in May in Fig. 2b. A larger vertical melt extent is therefore expected for this grid point. Moreover, the addition of internal heating enhances subsurface temperatures even more, further increasing the vertical melt extent. As a warmer snowpack is more prone to melting, more 305 melt events are expected to occur in the accumulation zone like in SE GRL. This is discussed in Sect. 6.1.
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 310 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 (Loewe, 1970). Here, we analyze whether this assumption is true given the observed impact of radiation penetration on summer subsur-315 face temperatures. Figure 5 shows the averaged yearly mean T 10m at the end of the simulation for Rp3 (Fig. 5a), the difference between the average skin temperature (T skin ) and T 10m (Fig. 5b) and the T skin − T 10m difference for Rp3 and Rp3 WIE with respect to Rp2 ( Fig. 5c and d, respectively).
In the interior, the difference between T skin and T 10m is small (< 0.5 • C) 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. 5c), for Rp3 320 WIE, however, it is considerably larger (often >2 • C, Fig. 5d). 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. 5d highlights the importance of radiation penetration.
In the percolation zone, the difference between T skin and T 10m is larger (more than -5 • C, Fig. 5b), owing to latent heat being released upon refreezing. Consequently, T 10m is not representative for the surface temperature in this zone. There are 325 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.

Surface energy budget
330 Figure 6 shows that the SEB components of Rp3 correlate generally well with observations of S6, S9 and S10 of the Ktransect and KAN-U and KAN-M of PROMICE. The downward and upward longwave radiative flux (LWd and LWu) show a determination coefficient (R 2 ) of 0.93 and 0.97, respectively, and the performance of Rp3 is similar to Rp2 (Table 1). The influence of cloud cover is captured well, with a small and insignificant bias of LWd and a small shortwave downward radiative flux (SWd) bias. The bias of SWd, however, is still significant, as is determined by statistical bootstrapping, due to the large 335 amount of data points. Despite its statistical significance, note that the uncertainty in the measurements (Sect. 2.5) is larger than the bias. 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 (green dot). The star indicates a location that is discussed in Sect. 4.1.
(d) The difference between Rp3 WIE and Rp2 for T skin − T10m. In (c) and (d), a positive value indicates that the T skin − T10m has become larger compared to Rp2. The upward shortwave radiative flux (SWu) is also in fairly good agreement with observations. The bias, however, is still statistically significant. The bias and RMSE of these fluxes have improved in Rp3 compared to Rp2 (Table 1), confirming that the albedo has improved, which is in agreement with Van Dalum et al. (2020).

340
Both the sensible and latent heat flux (SHF and LHF, respectively), show a large spread and significant bias ( 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 ) 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 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 longwave radiation (LWn), defined as LWd -LWu, is lower during winter for both Rp2 and Rp3, which is mostly due to the temperature bias induced by an overestimation of SHF. No changes have been made to the cloud and SHF parameterizations in Rp3, so similar results as Rp2 are therefore expected for cloudy conditions (Noël et al.,350 2018b). Improving the representation of surface roughness may also improve these results.
During summer, the SHF is overestimated by Rp2 and to a lesser degree by Rp3 for S6 (Fig. 7a). Refreezing, however, has increased for melting bare ice in Rp3 (discussed in more detail in Sect. 6.1 and Fig. 9), inducing more latent heat release, consequently heating the ice. As a result, the average skin temperature of the surface during the summer months increases and the temperature difference with the atmosphere reduces, leading to a smaller SHF. This effect is not present at S9 and S10, as 355 bare ice rarely reaches the surface at these sites.
During summer, the differences with observations increase for S6 and S9, as both model versions overestimate SWn due to a lower albedo compared to AWS observations. Note, however, that S6 is located in rough and inhomogeneous terrain, and that the local albedo and SEB may not be representative for the entire RACMO2 grid point.

360
Still, some differences between Rp2 and Rp3 are worth mentioning. At the onset of the accumulation season at S6, thin snow layers form and are modeled on top of bare ice. The lack of radiation penetration and internal heating in Rp2 leads to a too rapid brightening of the surface in Rp2, and subsequently underestimated SWn and hence melt-albedo feedbacks in early fall.
Eventually, the new snow layer becomes thick enough to cover the bare ice even if radiation penetration is taken into account, and the differences between both model versions disappear. This effect is also to a lesser degree visible at S9. For S9, which is 365 located close to the percolation zone, the albedo is overestimated for both Rp2 and Rp3. This is likely caused by an incorrect representation of the firn layer. The albedo differences are small for S10 (Van Dalum et al., 2020). As S10 is located in the lower accumulation zone, it is typically characterized by a homogeneous thick firn layer and little melt. Consequently, SWn is relatively similar in both model versions and is often in agreement with observations (Van Dalum et al., 2020).
6 Surface mass balance 370 6.1 Comparison with RACMO2.3p2 Figure 8a shows the annual average 2006-2015 SMB for Rp3. The absolute and relative difference with Rp2 (Rp3 -Rp2) is shown in Fig. 8b and c, respectively. Fig. 8b also shows the statistical significance. Almost all SMB changes are significant with respect to the inter-annual variability, only some of the smaller changes in the interior are insignificant. Integrated over the ice sheet, the SMB increases by 17.0 mm w.e. yr −1 (29.1 Gt yr −1 ), or 9.7%. The annual difference between Rp3 and 375 Rp2 and significance for melt, refreezing and runoff are shown in Fig. 9. Changes in modeled precipitation, sublimation and drifting snow erosion are not shown, as these processes are not significantly affected by the changes implemented in Rp3, and the differences are on average below 1.1 mm w.e. yr −1 . In the interior, a significant melt increase is observed in northeast and south Greenland (Fig. 9a). This is mostly caused by albedo changes, but also partly by radiation penetration inducing internal heating. As is discussed in Sect. 3, internal heating 380 raises the subsurface temperature and increases melt. Because all melt refreezes in the snow pack, these changes do not affect the local SMB ( Fig. 9b and c). Refreezing, however, changes the structure of the snow pack, inducing more energy absorption and higher snow temperatures, leading to more melt.
The outer rim of the ice sheet, except in the southeast, is characterized by a strong and significant SMB increase (Fig 8).
In Rp2, the ice mask is not projected properly on the bare ice albedo field, causing the outermost glaciated grid points to be 385 contaminated with tundra albedo. The albedo of the outermost glaciated grid points is therefore as low as 0.3, inducing too much melt and runoff. The impact of this artifact is mitigated for higher resolutions, and is solved in Rp3 (Van Dalum et al., 2020), lowering the 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. Melt is slightly higher and refreezing is similar to Rp2 for this region, resulting in more runoff.

390
The lower ablation area of the GrIS in the southwest experiences a strong increase in refreezing (Fig. 9b), while the SMB difference is limited. Refreezing is enhanced by the introduction of radiation penetration. In Rp3, subsurface melting of ice creates pore space (Sect. 2.1.1) and consequently increases the melt water retention capacity. In Rp2, no retention capacity remains once bare ice reaches the surface, limiting refreezing in the shallow winter snowpack. Additionally, some regions in the southwest show a melt increase, which is induced by, on average, a slightly lower albedo in Rp3 (Van Dalum et al., . As a result, the increase in melt and refreezing balance out, leading to an insignificant runoff difference for most of the southwestern ablation zone (Fig. 9), hence the SMB changes little (Fig. 8b). and cloudy conditions. A minimum vertically integrated cloud content of 0.05 kg m −2 is required to be considered cloudy; and a minimum melt rate of 1 mm w.e. day −1 is required to be considered a melt event. Using this low threshold should guarantee that all major melt events are included (Fettweis et al., 2011;Fausto et al., 2016). All data are averaged between 2006 and 2015.
In Rp3, 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. 8c). Similarly, in the upper ablation zone, where the winter snow cover lasts long during 400 summer, enhanced melt leads to reduced refreezing, as the refreezing capacity is consumed faster.
The northeast shows a large area with reduced SMB (pattern B in Fig. 8c). 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 prolonging the melt season (as is discussed in Sect. 4.1), consequently enhancing melt and runoff (Fig. 9a, c, d and f).

405
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.
The larger melt extent for the star in south-east Greenland (Fig. 4c, star in Fig. 5c) is in agreement with Fig. 9d and e, where we see more melt and refreezing for this grid point. As the elevation of this location is too far above sea level, no runoff is 410 modeled and all melt water refreezes locally. Even though the melt and consequent refreezing are still small and do not alter the SMB, they do change the snow structure. for a large area in the center (Fig. 9d and e). This is discussed in more detail in Sect. 6.2. No runoff is modeled in the interior ( Fig. 9c and f), as melt water refreezes locally. Around the K-transect, the increase of refreezing consequently lowers the modeled runoff. Significant runoff increase is modeled in A and B.

Cloud cover and melt 420
For a large area of central GrIS, annual melt and refreezing have decreased by percentage ( Fig. 9d and e). Melt events in this region, although rare, are almost always associated with above-average cloud cover (Fig. 10a), 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., 2020), but for cloudy conditions and during melt events, the albedo of Rp3 is higher (Fig. 10b), 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 425 distribution of irradiance during cloudy conditions, as relatively more of the 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 (Dang et al., 2015;Van Dalum et al., 2019;Warren, 2019). A large fraction of UV and visible radiation is absorbed internally (Fig. 3), heating subsurface layers that can potentially increase melt. This effect, however, is not enough to mitigate the reduced irradiance and higher albedo, and therefore less melt is modeled. To summarize, cloud cover during melt events changes the spectral distribution of shortwave 430 radiation at the surface, shifting more towards shorter wavelengths (UV and visible light). As the albedo is higher for these shorter wavelengths, which is now properly taken modeled in Rp3, less energy is available for melt. Note that the total amount of 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).

SMB Observations
435 Figure 11 shows the annual SMB of Rp3 and Rp2 with respect to various SMB observations. The determination coefficient, root-mean-square error (RMSE) and bias with respect to observations are also shown. A slightly higher SMB is modeled for Rp3 compared to Rp2. The bias with observations for Rp3 is -0.031 m w.e. yr −1 (-2.67%) and for Rp2 -0.091 m w.e. yr −1 (-7.54%), but both are statistically insignificant. The differences with observations are particularly relevant for the measurement sites in the lower ablation zone (S6, squares in Fig. 11) and are caused by increased refreezing in Rp3 (Fig. 9b). Note that 440 the spread is 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 close to the ice margin.

SMB without internal energy absorption
Finally, we discuss the impact of radiation penetration on the SMB and analyze its components by comparing Rp3 with Rp3 445 WIE. Figure 12 shows the differences in the annual SMB, melt and refreezing between the two simulations. 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 is higher for Rp3 in most of the ablation zone (Fig. 12a), in particular in the southwest (area C of Fig. 12c), where we observe a relatively small melt decrease (Fig. 12b). More importantly, however, is that internal radiative heating in Rp3 increases refreezing in the ablation zone (Fig. 12c), especially in the southwest, but also in the northwest and northeast (area 450 C, D and B of Fig. 8c, respectively), reducing modeled runoff. Hence, the SMB of Rp3 WIE is lower than Rp3 for these areas.
This increase in refreezing in Rp3 is likely realistic, as bare ice does have retention capacity for water and refreezing can occur overnight, while in Rp2 melting bare ice essentially remains dry as meltwater is removed instantaneously.
In the central part of the interior, the melt and refreezing differences between Rp3 and Rp3 WIE are smaller than for the surrounding areas (Fig. 12b, c). This supports the findings of Sect. 6.2, where we show that the albedo dominates the signal 455 and not subsurface energy absorption. The differences between Rp3 and Rp3 WIE are consequently only small, as both model simulations use the same albedo scheme. For the surrounding areas, however, Fig. 12b and c shows that subsurface energy absorption plays a more important role, as the difference between the model simulations for melt and refreezing increases.
In south Greenland, the melt and refreezing difference between Rp3 and Rp3 WIE is smaller than between Rp3 and Rp2 ( Fig. 9d). This means that processes other than subsurface energy absorption are important here. This region is characterized 460 by an albedo decrease due to slightly stronger snow metamorphism and slower firn-ice transition in Rp3 (Van Dalum et al., 2020), resulting in more absorption of shortwave radiation, inducing melt and subsequent refreezing. Hence, it contributes to the differences between Rp3 and Rp2.

Conclusions
Regional climate models often do not explicitly account for radiative transfer in snow and ice. Therefore, we evaluated the 465 latest version of the regional climate model RACMO2, which has an updated radiative transfer scheme that allows internal heating. In this new version Rp3, we also updated the multilayer firn module and snow and ice albedo (Van Dalum et al., 2020), and developed a method to distribute total absorbed shortwave radiation between the SEB and to internal absorption. In this study, we assessed the modeled surface mass balance (SMB), surface energy budget (SEB) and firn and snow temperature of the GrIS.
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 temperature profiles show a more gradual temperature decrease with depth in Rp3 during summer in the interior, improving the agreement with observed snow temperatures at Summit; and melt reaches deeper layers both in the accumulation and ablation zone.

475
Compared to the previous RACMO2 version Rp2, the SMB difference is small in the interior, as any change in melt is balanced by refreezing. Close to the percolation zone near 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, significantly 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 melt. These differences are 480 small or absent in the wet southeast, where bare ice rarely reaches the surface. For the accumulation zone of the southeast, relative differences in melt and refreezing are substantial due to subsurface energy absorption and melt water percolation. This increased melt and refreezing often does not lead to significant SMB changes, but does change the snow structure.
There is, however, still room for improvement. The resolution used in this study, i.e., 11 km horizontal resolution, is sufficient to resolve the SMB and SEB in the interior of the ice sheet, but is insufficient to resolve the inhomogeneous ablation zone, 485 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 horizontal resolution of RACMO2. More fundamentally, the ice albedo, although considerably improved in Rp3 (Van Dalum et al., 2020), could be further improved 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 space and time and would be beneficial. The turbulent fluxes could be further improved as well. 490 We have shown that there are still some biases in the SHF and LHF, especially for areas with rough terrain like S6. A better representation of the surface roughness and the atmospheric fluxes in this complex terrain is therefore desirable (Van Tiggelen et al., 2021).
To conclude, Rp3 is capable to accurately reproduce the SMB, SEB and (sub)surface temperature when compared to in-situ observations. Differences with Rp2 are most pronounced in the ablation zone close to the ice margin, in south Greenland and 495 close to the percolation zone. Furthermore, the snow temperature with depth is altered ice-sheet wide due to subsurface energy absorption. The improvements made in Rp3 are especially relevant in a warming climate where the GrIS surface will typically have lower albedos and more subsurface energy absorption, thus improving RACMO2's capability to make future climate projections. This study therefore shows the necessity of radiative transfer in snow and ice for regional climate modeling of the GrIS. In a subsequent study, we will assess the ability of Rp3 to simulate albedo and melt in Antarctica.
Rp2 data are available from the authors.
Competing interests. The authors declare that they have no conflict of interest.