3.1. Surface energy-balance model

The surface energy balance is calculated by:

(1)$$S_{\rm in} + S_{\rm ref} + L_{\rm in} + L_{\rm out} + H + LE + Q_{\rm R} + Q_{\rm G} + Q_{\rm M} = 0,\; $$

where $S_{\rm in}$ and $S_{\rm ref}$ are the incoming and reflected shortwave radiation, respectively, $L_{\rm in}$ and $L_{\rm out}$ are the incoming and outgoing longwave radiation, respectively, $H$ is the sensible heat flux, $LE$ is the latent heat flux, $Q_{\rm R}$ is the energy supplied by rain, $Q_{\rm G}$ is the energy flux into the subsurface (also referred to as ground heat flux) and $Q_{\rm M}$ is the energy available for melt. Fluxes toward the surface are defined as positive. The model uses a skin layer formulation to close the energy balance at the surface. In this approach the surface is assumed to be an infinitesimal skin layer without heat capacity that reacts instantaneously to a change in energy input. Several fluxes in Eqn (1) directly depend on the surface temperature $T_{\rm s}$ ($L_{\rm out}$, $H$, $LE$, $Q_{\rm G}$ and $Q_{\rm R}$). A bisection method is iteratively applied to find the $T_{\rm s}$ that satisfy Eqn (1) assuming $Q_{\rm M} = 0$. If a positive surface temperature is found (indicating surface melt, i.e. $Q_{\rm M} \ne 0$), $T_{\rm s}$ is reset to zero, and the energy balance is recalculated. The resulting non-zero $Q_{\rm M}$ is converted into melt.

The radiative fluxes $S_{\rm in}$, $S_{\rm ref}$ and $L_{\rm in}$ are taken from observations at the weather stations. $L_{\rm out}$ is calculated from the modeled surface temperature using the Stefan–Boltzmann law, assuming an emissivity of $1$. Measured $L_{\rm out}$ is used for model validation.

The turbulent fluxes are calculated using the bulk aerodynamic method, relating $H$ and $LE$ to wind speed and gradients of temperature and vapor pressure between the surface and the air (Hock and Holmgren, Reference Hock and Holmgren2005). Wind speed and air temperature are based on observations, and vapor pressure is calculated from measurements of relative humidity. For these calculations, we assume a constant measurement height above the surface of 2 m. The effect of atmospheric stability is taken into account using the Monin–Obukhov similarity theory. Stability functions for stable atmospheric stratification are computed following Beljaars and Holtslag (Reference Beljaars and Holtslag1991), while stability functions for unstable atmospheric stratification are computed following Panofsky and Dutton (Reference Panofsky and Dutton1984). Surface roughness length for momentum is set at a constant value of $10^{-4}$ m, which is found to be a good assumption at these elevations on the ice sheet (Smeets and Broeke, Reference Smeets and Broeke2008a; Charalampidis and others, Reference Charalampidis2015). Surface roughnesses lengths for heat and moisture are calculated using the surface renewal theory (Andreas, Reference Andreas1987; Munro, Reference Munro1990).

The energy supplied by rain $Q_{\rm R}$ is determined by the rainfall rate and the temperature difference between the air and the surface (Hock and Holmgren, Reference Hock and Holmgren2005). The ground heat flux $Q_{\rm G}$ is calculated from the temperature gradient between the surface and the first subsurface layer using Fourier's law and the thermal conductivity provided by the subsurface model. In this calculation, modeled surface temperature from the current time step and subsurface temperature from the previous time step are used.

3.2. Subsurface model

The subsurface model is based on the SOMARS model developed by Greuell and Konzelmann (Reference Greuell and Konzelmann1994) and further refined by Reijmer and Hock (Reference Reijmer and Hock2008). The model calculates firn temperature, density of the dry part of the firn and liquid water content on a vertical grid, including the effect of meltwater percolation, refreezing and dry firn densification. Here the grid extends from the surface to a depth of $25$ m with a total of $\sim$$140$ layers with varying resolution: $0.05$ m from the surface to $1$ m depth, $0.1$ m between $1$ and $10$ m depth and $0.5$ m below $10$ m depth.

First, the temporal evolution of the firn temperature is calculated using a forward time centered space finite difference method to solve the 1-D heat equation. Boundary conditions are given by the surface temperature at the snow surface, computed by the surface energy balance model, and by a constant temperature at the bottom of the firn column, determined by observations ($-13.4^{\circ }$C at Site J and $-19.3^{\circ }$C at EKT). The thermal conductivity ($\kappa$) is described as a function of the firn density, following the approach presented in Douville and others (Reference Douville, Royer and Mahfouf1995).

Next, if energy is available for melting ($Q_{\rm M} > 0$), the amount of melt at the surface is calculated. Meltwater and rainwater are then added to the liquid water content of the first subsurface layer. Refreezing, percolation and densification are sequentially calculated for every layer. The amount of refreezing is limited by either layer temperature, liquid water content or available pore space. If any water is still present in the layer, a small amount of it is retained in the layer against gravity due to capillary forces (irreducible water content) while the rest is allowed to percolate to the next layer. The irreducible water content is described as a function of the firn density, according to Schneider and Jansson (Reference Schneider and Jansson2004). Finally, the densification of the dry firn is calculated depending on the temperature gradient and the density following the Herron and Langway (Reference Herron and Langway1980) parameterization adapted by Li and Zwally (Reference Li and Zwally2004).

Solid precipitation is added to the first layer using observations of fresh snowfall thickness and density. The thickness for each time step was derived from total daily surface height changes in order to filter the often noisy sonic ranger data; daily values were then distributed equally over the day for the simulations. A fresh snowfall density of $350$ kg m$^{-3}$ was derived from snow pit observations at Site J and EKT in Spring $2017$. This value reflects the rapid, wind-driven, compaction following snowfall events and it is in line with previous studies suggesting values between $315$ and $400$ kg m$^{-3}$ for this area of the ice sheet (Charalampidis and others, Reference Charalampidis2015; Fausto and others, Reference Fausto2018; Vandecrux and others, Reference Vandecrux, Fausto, Langen, van As and Macferrin2018). The model does not account for the lateral flow of water. Subsurface layers are allowed to merge or split to maintain the initial grid configuration if they become thinner than $50$% or thicker than $150$% of their original thickness. This is often the case for the first layer, where thinning and thickening occur during melting and snowfalls.

3.3. Model application

The model was run at Site J and EKT from $6$ May $2017$ to $4$ September $2019$ for all the periods when meteorological data were available (Fig. 1c). Hourly data of air temperature, relative humidity, wind speed, shortwave incoming and reflected radiation, longwave incoming radiation and solid precipitation provided the meteorological input to run the model. All the forcing inputs were taken directly from weather station observations. The model was run at $4$ min resolution to assure the stability of the subsurface solver; hourly input data were linearly interpolated to match the required time step. The subsurface model was initialized with firn temperature and density data from the firn thermistor string and firn cores, respectively. When firn core data were unavailable (May $2018$ at Site J and EKT, May $2019$ at Site J), the firn density profile was initialized using model output from the end of that site's previous simulation adjusted to account for the correct seasonal snow depth.

3.4. Model validation

Simulations are validated against independent observations of relative surface height from both the weather station sonic ranger and ablation stakes readings (Fig. 2), surface temperature estimated from the measurements of outgoing longwave radiation using the Stefan–Boltzmann law with an emissivity of 1 (Fig. 2), and subsurface temperatures from thermistor strings (Fig. 3).

Fig. 2. Model validation at EKT (a, b) and Site J (c, d). (a, c) Modeled hourly relative surface height compared to measurements from the AWS sonic ranger (SR50) and ablation stake between May $2017$ and September $2019$. (b, d) Modeled vs observed hourly surface temperature, with the 1:1 line in black. ${N}$ is the number of samples, RMSE the root mean square error, ${r}$ the correlation coefficient and ${p}$ the $p$-value.

Fig. 3. Hourly (a, d) measured, (b, e) modeled subsurface temperatures and (c, f) their differences for the uppermost 10 m between May $2017$ and September $2019$ at EKT (a–c) and Site J (d–f). The depths are relative to the snow surface on the date of installation in early May $2017$. Black lines and black dots in (a, c, d, f) indicate the observed relative surface height and the thermistors position, respectively. Blue lines in (b, e) contour the boundary of subsurface liquid water. Gray lines in (b, e) indicate the depth of the uppermost thermistor. Data from the new thermistors installed in May $2019$ are omitted in (a, c, d, f) when the sensors reached the surface or were affected by solar radiation penetration.

Independent in situ ablation stake readings, which provide relative surface height changes between annual field visits, match the continuous sonic ranger (SR50) data very well (Figs 2a, c). This provides confidence that the weather stations, where the sonic rangers were mounted, were not subject to sinking during the study period. In contrast, simulated relative surface heights do not agree well with these observations. Modeled surface lowering, almost entirely caused by melt, is consistently overestimated at both sites and during all three years except for EKT in summer 2017, where only little melt occurred. Maximum differences between observed and simulated relative surface height up to 0.47 m at Site J and 0.42 m at EKT are recorded during summer $2019$ when surface melting was most pronounced.

The comparison of modeled hourly surface temperatures with those derived from the measurements of outgoing longwave radiation reveals that the model tends to overestimate $T_{\rm s}$, especially during the summer months when melt occurs at the surface (Figs 2b, d). The average difference between modeled and observed temperatures is $+ 1.11^{\circ }$C at EKT and $+ 0.59^{\circ }$C at Site J over the whole study period, increasing to $+ 1.93$ and $+ 1.31^{\circ }$C, respectively when only the summer months (June, July and August) are considered. The correlation coefficient between hourly model results and observations is $0.97$ and $0.99$, and root mean square error is $2.8$ and $1.9^{\circ }$C at EKT and Site J, respectively.

Figure 3 shows simulated and observed subsurface temperatures and their differences. Temperatures below $-10$ m depth (not shown) remain constant at both sites and model results agree well with the observations. Simulated temperatures generally match the observations at depths below $4$ m well, but the model tends to underestimate near-surface (between $0$ and $4$ m depths) temperatures during the winter and overestimate them during the summer (Figs 3c, f). The overestimation of summer temperatures is more evident in $2019$ when melt was more pronounced. Furthermore, in May $2019$, two new thermistors were installed at both sites in that year's winter snow ($\sim$first meter), allowing for better validation of near-surface model results. Such observations were not available in $2017$ and $2018$ (first thermistors at $-1$ m or lower), preventing us from validating model results closer to the surface in those years.

At both sites, modeled meltwater percolation (contoured by blue lines) did not reach the depth of the uppermost thermistor (gray lines) during the $2017$ and $2018$ melt seasons (Figs 3b, e). During summer $2019$, modeled meltwater percolated much deeper than in the previous summers, reaching $-1$ m at EKT and $-4$ m at Site J. However, observations do not show the same thermal signal of percolating meltwater, typically characterized by the snowpack and firn becoming isothermal ($0^{\circ }$C). This results in large temperature differences (Figs 3c, f), with modeled temperatures up to $+ 13.9$and $+ 10.8^{\circ }$C warmer at EKT and Site J, respectively, during summer $2019$.

4.1. Model forcings

The meteorological observations used to force the model can be prone to errors. Here we test how such errors could affect model results by homogeneously perturbing the five meteorological input variables one by one, and re-running the model with the new input forcings (Fig. 4). Applied perturbations are guided by instrument accuracy as declared by the manufacturer (Table I). For example, we perturb air temperature by $\pm 1^{\circ }$C when the manufacturer stated accuracy is $0.6^{\circ }$C, etc. Furthermore, for each of the input variables, we explore a scenario that, on average, comes close to relative surface height observations (shown in orange in Fig. 4).

Fig. 4. Modeled hourly relative surface height compared to measurements from AWS sonic ranger (SR50) and ablation stake at Site J for simulations with different model forcings: (a) air temperature, (b) relative humidity, (c) wind speed, (d) net shortwave radiation and (e) incoming longwave radiation.

These experiments suggest that the discrepancies between the reference simulation and observations cannot be explained by errors in input forcing or instrument accuracy. The perturbations required to match the measured surface height are extreme and far from being physically reasonable, e.g.$-5^{\circ }$C to air temperature, $-50$% to relative humidity, etc. Finally, given the large differences between model results and measurements, we are confident that not even the effect of combined perturbations (e.g. perturbing simultaneously two or more variables within reasonable limits) could explain the discrepancies.

Figure 4 indicates that, while positive perturbations typically produce more surface lowering, this is not the case for the wind speed. In fact, because turbulent fluxes at our sites are generally negative, higher wind speeds would cause more negative fluxes resulting in less surface lowering.

4.2. Model parameters and parameterizations

Model parameters and parameterizations (e.g. turbulent fluxes, thermal conductivity, etc.) used in the reference simulation are taken from the literature rather than calibrated to match the observations. To further investigate possible causes for the mismatch between model results and observations in $sim\_ref$, here we test the model sensitivity to different choices of model parameters and parameterizations.

We perturb the surface roughness for momentum used in the reference simulation for the calculation of the turbulent heat fluxes ($z_{0m} = 0.1$ mm) by increasing and decreasing it by one order of magnitude (Fig. 5a). We change the new snow density from $350$ kg m$^{-3}$ to $300$ and $400$ kg m$^{-3}$ (Fig. 5b). Thermal conductivity ($\kappa$) determines subsurface temperature evolution; thus, it controls the ground heat flux calculation, ultimately affecting the surface energy balance. Several different parameterizations for thermal conductivity exist in the literature. Here we replace the one used in the reference simulation (Douville and others, Reference Douville, Royer and Mahfouf1995) with four others following Reijmer and Hock (Reference Reijmer and Hock2008): Van Dusen and Washburn (Reference Van Dusen and Washburn1929); Sturm and others (Reference Sturm, Holmgren, König and Morris1997); Östin and Andersson (Reference Östin and Andersson1991); Jansson (Reference Jansson1901) (Fig. 5c). Similarly, we replace the irreducible water content parameterization by Schneider and Jansson (Reference Schneider and Jansson2004) used in the reference simulation with the one presented in Coléou and Lesaffre (Reference Coléou and Lesaffre1998) (Fig. 5d). Finally, we turn off the densification parameterization (Fig. 5d).

Fig. 5. Model sensitivity to choice of model parameters and parameterizations at Site J. Modeled hourly relative surface height compared to measurements from sonic ranger (SR50) and ablation stake for simulations with (a) perturbed roughness length for momentum and (b) perturbed fresh snow density, (c) different thermal conductivity and (d) densification and irreducible water content ($\Theta _{mi}$) parameterizations.

As evident from Figure 5, the model is not sensitive to these perturbations and changes. Differences to $sim\_ref$ are barely detectable in most cases. Changes in surface roughness for momentum by an order of magnitude only slightly change the simulations. It should be noted that because the model uses the parameterization by Andreas (Reference Andreas1987), surface roughness lengths for heat and moisture are computed as a function of $z_{0m}$; thus, they are also affected by these perturbations.

While not directly affecting the surface energy balance, the fresh snow density parameter controls the conversion of melt into surface height change. Higher densities cause less surface lowering and vice versa. This is only noticeable when significant accumulation occurred before or during the melt season, e.g. in $2017$ and $2018$ at Site J (Fig. 5b). Overall, simulation results are barely affected by the experiments. Thus, uncertainties in roughness lengths, snow density and thermal conductivity cannot explain the discrepancies between the reference simulation and observations.

4.3. Penetration of shortwave radiation

Although a well-documented physical process (e.g. Schlatter, Reference Schlatter1972; Colbeck, Reference Colbeck1989; Brandt and Warren, Reference Brandt and Warren1993), the penetration of shortwave radiation into the snowpack is not included in our model, consistently with most state of the art surface energy-balance models, including the most recent versions of RCMs such as MARv3.12, HIRHAM5 and RACMOv2.3p2 (Mankoff and others, Reference Mankoff2021). Few studies have investigated the effect of shortwave radiation penetration on the surface energy balance in Greenland (Van Den Broeke, Reference Van Den Broeke2008; Kuipers Munneke and others, Reference Kuipers Munneke2009) and in the Antarctic Peninsula (Kuipers Munneke and others, Reference Kuipers Munneke, van den Broeke, King, Gray and Reijmer2012). Here we use the same approach as Kuipers Munneke and others (Reference Kuipers Munneke2009), described in detail in the Appendix, to assess whether including this process in our model could explain the differences from the observations presented above.

The penetration of shortwave radiation, as described in this parameterization, strongly depends on the thickness of the fictitious surface layer, which is not to be confused with the thickness of the subsurface layers (see the Appendix for details and Kuipers Munneke and others (Reference Kuipers Munneke2009)). This calibration parameter controls the amount of net shortwave radiation ($S_{\rm net}$) that penetrates into the subsurface ($S_{\rm pen}$) and the amount that is absorbed by the surface ($S_{\rm sfc}$), affecting the surface energy balance and the warming of the subsurface. Here we present simulations using three different fictitious surface layer thicknesses: $0.01$, $0.05$ and $0.25$ m (Fig. 6). For small thicknesses, more shortwave radiation is allowed to penetrate into the subsurface while for thicker layers most of the shortwave radiation is absorbed by the surface: $S_{\rm pen} / S_{\rm net} = 43$, $17$ and $4$% on average for the $0.01$, $0.05$ and $0.25$ m simulations respectively (Fig. 6c). If a very thick fictitious surface layer is considered all the net shortwave radiation would stay at the surface, replicating the assumptions made in the reference simulation when radiation penetration is not considered.

Fig. 6. (a) Hourly relative surface height, (b) annual cumulative total, surface and subsurface melt and (c) $30$-day rolling average of shortwave radiation components ($S_{\rm net}$, $S_{\rm sfc}$ and $S_{\rm pen}$) for simulations including radiation penetration at Site J. Results using different fictitious surface layer thicknesses are presented: $1$ cm in cyan, $5$ cm in orange and $25$ cm in red. Reference simulation in solid black and observations in dashed black. Reference simulation only has total cumulative melt and $S_{\rm net}$.

Both relative surface height and total melt vary very little from $sim\_ref$ (Figs 6a, b). This is because radiation penetration allows melting to happen also in the subsurface (e.g. firn temperature $= 0^{\circ }$C). Hence, the thinner the fictitious surface layers, the greater the melt in the subsurface layers, thereby compensating for less melt at the surface (Fig. 6b).

More details on the effects of radiation penetration are highlighted when looking at surface and subsurface temperatures. Figure 7 shows hourly temperatures for a 2 weeks period at Site J during the $2019$ melt season. During the first days, the near-surface snowpack is mostly temperate (i.e. at $0^{\circ }$C), indicating the presence of meltwater. Starting $28$ June the air temperature drops resulting in a cooling of both the surface and the subsurface and an almost complete halt in melt. Surface temperature is higher for simulations with thicker fictitious surface layers, and it is the highest for the reference simulation (Fig. 7a), where radiation penetration is not considered. However, the simulated $T_{\rm s}$ is higher than the observations for all the simulations. In contrast, subsurface temperatures, at both $0.1$ and $0.25$ m depth, are warmer for simulations with a thinner fictitious surface layer when more energy penetrates into the subsurface. Despite subsurface temperatures being sensitive to the choice of the fictitious surface layer thickness necessary to model the penetration of solar radiation, the experiments indicate that the omission of this process in the reference simulation cannot account for the differences between model results and observations.

Fig. 7. Details of hourly surface temperature (a) and subsurface temperature at $0.10$ m (b) and $0.25$ m (c) depth for simulations including radiation penetration at Site J for a 2 weeks period between the end of June and the beginning of July $2019$. Results using different fictitious surface layer thicknesses are presented: $1$ cm in cyan, $5$ cm in orange and $25$ cm in red. Reference simulation in solid black and observations in dashed black.

4.4. Deep water percolation

Deep water percolation through preferential flow paths is a well-known process that affects water redistribution through the snow and firn matrix (e.g. Sturm and Holmgren, Reference Sturm and Holmgren1993; Humphrey and others, Reference Humphrey, Harper and Pfeffer2012; Cox and others, Reference Cox, Humphrey and Harper2015). However, this process is often not included in the subsurface models used in surface energy-balance models, which employ a tipping bucket approach to simulate water percolation from one layer to the next one. This is also the case for RCMs such as MAR, HiRAM and RACMO2. Marchenko and others (Reference Marchenko2017) used a simple statistical approach to parameterize deep water percolation, showing that accounting for this process improves the simulations of subsurface temperature.

Here, we apply the same approach, described in detail in the Appendix, to assess whether parameterizing deep water percolation could improve our model results presented above. Surface meltwater is redistributed to deeper subsurface layers following a probability density function that solely depends on a prescribed depth of maximum percolation (e.g. depth below which no water is allocated). Following Marchenko and others (Reference Marchenko2017), we test three different probability density functions, uniform (UNI), linear (LIN) and normal law (NORM), and three different maximum percolation depths, $2.5$, $5.0$ and $7.5$ m.

Accounting for deep water percolation reduces the $sim\_ref$ total melt at Site J over the entire study period between $1.4$ and $9.7$%, depending on the probability density function and the maximum percolation depth used (Table S1). Surface melt is reduced the most when more water is redistributed from the surface to deeper subsurface layers (Fig. 8b). This is the case for simulations with deeper maximum percolation depths (e.g. $7.5$ m) and using functions that allocate more water at depth (e.g.UNI). However, despite the decrease in melt, the relative surface height lowers more in these simulations than in $sim\_ref$ (Fig. 8a). This is the result of redistributing meltwater and refreezing away from the surface. In a tipping bucket percolation model, in fact, the refreezing capacity of a layer is completely exhausted before water is allowed to percolate to the next layer. Thus, in $sim\_ref$ refreezing rapidly increases the density of the surface layer making the surface lowering less compared to simulations with deep percolation. On the contrary, when water is allowed to percolate at depth, refreezing happens at multiple layers delaying the appearance of liquid water in the subsurface (Fig. S1) and keeping the density of the surface layer low.

Fig. 8. (a) Hourly relative surface height and (b) annual cumulative melt for simulations including deep water percolation at Site J. Results using different percolation probability density functions (UNI in cyan, LIN in orange and NORM in red) and depth of maximum percolation ($2.5$ m dotted, $5.0$ m solid and $7.5$ m dashed) are presented. Reference simulation in solid black and observations in dashed black.

The redistribution of refreezing at depth that follows deep water percolation also affect the subsurface temperature. All simulations overestimate the subsurface temperature at depths where water is allowed to percolate, showing larger temperature differences between model and observations compared to $sim\_ref$ (Fig. S2). Temperature differences are the largest for simulations that redistribute more water at depth, which are the same that also reduce melt the most.

These experiments indicate that, although accounting for deep water percolation slightly decreases melt, it makes the comparison with observations even worse than for the reference simulation.