Modelling snowpack on ice surfaces with the ORCHIDEE land surface model: application to the Greenland ice sheet

Current climate warming is accelerating mass loss from glaciers and ice sheets. In Greenland, the rates of mass changes are now dominated by changes in surface mass balance (SMB) due to increased surface melting. To improve the future sea-level rise projections, it is therefore critical to have an accurate estimate of the SMB, which depends on the representation of the processes occurring within the snowpack. The Explicit Snow (ES) scheme implemented in the land surface model Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) has not yet been adapted to ice-covered areas. Here, we present the preliminary developments we made to apply the ES model to glaciers and ice sheets. Our analysis mainly concerns the model's ability to represent ablation-related processes. At the regional scale, our results are compared to the MAR regional atmospheric model outputs and to MODIS albedo retrievals.

Using different albedo parameterizations, we performed offline ES simulations forced by the MAR model over the 2000–2019 period. Our results reveal a strong sensitivity of the modelled SMB components to the albedo parameterization. Results inferred with albedo parameters obtained using a manual tuning approach present very good agreement with the MAR outputs. Conversely, with the albedo parameterization used in the standard ORCHIDEE version, runoff and sublimation were underestimated. We also tested parameters found in a previous data assimilation experiment, calibrating the ablation processes using MODIS snow albedo. While these parameters greatly improve the modelled albedo over the entire ice sheet, they degrade the other model outputs compared to those obtained with the manually tuned approach. This is likely due to the model overfitting to the calibration albedo dataset without any constraint applied to the other processes controlling the state of the snowpack. This underlines the need to perform a “multi-objective” optimization using auxiliary observations related to internal snowpack processes. Although there is still room for further improvements, the developments reported in the present study constitute an important advance in assessing the Greenland SMB with possible extension to mountain glaciers or the Antarctic ice sheet.

1 Introduction

Satellite observations reveal that the Greenland ice sheet (GrIS) has been losing mass for at least 3 decades. Between 1992 and 2018, the net ice mass loss was estimated at 3800±339 Gt, corresponding to a rise in global mean sea level of 10.6±0.9 mm (The IMBIE team, 2020). Mass loss is driven by dynamic solid ice discharges (Enderlin et al., 2014) and by enhanced surface meltwater and runoff (Ryan et al., 2019). Over the 2000–2008 period, the GrIS mass loss was equally partitioned between surface and dynamic processes (van den Broeke et al., 2009). However, recent studies based on regional climate models and remote-sensing observations (van den Broeke et al., 2016; Ryan et al., 2019; The IMBIE Team, 2020; Fox-Kemper et al., 2021) show that rates of mass change are now dominated by changes in surface mass balance (SMB), defined as the difference between mass gains (solid and liquid precipitation) and surface ablation processes (runoff, sublimation, and snow erosion).

Besides directly impacting the global mean sea level, the GrIS is also an integral part of the Earth system (Fyke et al., 2018). As such, it is highly sensitive to climate change and in turn has a strong influence on global climate, notably by releasing fresh water into the ocean, which leads to changes in the Atlantic meridional overturning circulation (Bakker et al., 2016; Martin et al., 2022). Surface melting may also induce changes in the local climate through temperature–elevation feedback (Edwards et al., 2014; Sellevod et al., 2019) and the albedo effect (Box et al., 2012; Helsen et al., 2017; Riihelä et al., 2019). Finally, changes in topography produce modifications in the local and large-scale atmospheric circulations (Ridley et al., 2005; Hahn et al., 2020).

To capture this feedback and to reduce the uncertainties in sea-level and climate projections, a key objective of the climate–ice-sheet modelling community is to incorporate ice-sheet models in Earth system models (ESMs; Vizcaino, 2014). Such coupled climate–ice-sheet models have mainly been developed with low-resolution climate models designed for long-term integrations (Kageyama et al., 2004; Charbit et al., 2005; Vizcaino et al., 2010; Roche et al., 2014). So far, only a few groups have met this goal with CMIP-like models (Vizcaino et al., 2013; Muntjewerf et al., 2020; Smith et al., 2021). A key challenge in developing such models relates to the realistic computation of SMB used as a forcing field of the ice-sheet models.

SMB is highly dependent on the radiative properties of snow and on the physical processes occurring within the snowpack (Helsen et al., 2017). At the surface, snow cover evolves as a function of the surface energy balance and mass exchanges with the atmosphere. In cold regions, snowmelt is largely driven by shortwave radiation: because of the high albedo value of fresh snow (0.80–0.90), a large fraction of shortwave radiation is reflected to the atmosphere, limiting the energy available at the surface for melting. Therefore, snow evolution is strongly dependent on the albedo. The value of snow albedo decreases when snow is ageing (i.e. in the absence of a new snowfall event) and with the snow metamorphism and liquid water content at the ice-sheet surface coming from either rainfall or snow/ice melting. Surface water may also percolate and refreeze inside the snowpack, thereby delaying the runoff. The transformation of snow into ice depends on environmental conditions (e.g. winds, near-surface temperatures) and internal processes within the snowpack (e.g. heat conduction and vertical temperature gradient, compaction), which directly influence the grain microstructure and the snow density. All these processes affect the SMB of the ice sheet.

There are several ways to compute the SMB. Empirical approaches such as the positive-degree-day method (Reeh, 1991) have long been used to compute snow and ice melting from downscaled near-surface temperatures. This kind of approach requires little in terms of computational resources and has often been applied to past and future long-term integrations (Charbit et al., 2008, 2013; Bonelli et al., 2009; Vizcaino et al., 2010). However, such methods have been calibrated against the present state of the GrIS, raising the question as to whether they can be applied to a different climatic context from the present-day one, knowing that ablation is projected to increase (van de Wal, 1996; Bougamont et al., 2007). Moreover, they are not physically based and cannot reproduce the diversity of snow processes that directly influence the SMB. Snow models implemented in general circulation models have long been based on simplified physics. They are mainly designed to resolve the seasonal and diurnal variations in heat fluxes but with no representation of internal processes (Armstrong and Brun, 2008). In contrast, regional climate models developed for polar regions generally incorporate multiple-layer energy balance snow models with a fine vertical resolution (e.g. Brun et al., 1992; Lefebre et al., 2003; Vionnet et al., 2012; Noël et al., 2018) and with detailed snow physics to simulate a variety of snowpack processes. However, due to their high computational cost, they are not often used in ESMs, despite a few rare exceptions such as the work of Punge et al. (2012) based on the implementation of a detailed snow model (Brun et al., 1992) in the atmospheric model LMDZ4 (Hourdin et al., 2006) or the Community Land Model (CLM), which includes the snow radiative-transfer scheme SNICAR (Flanner and Zender, 2006) and a snow model simulating a variety of key snow processes such as the metamorphism (Lawrence et al., 2019; He et al., 2024).

An alternative approach consists of implementing snow models of intermediate complexity in the land surface components of ESMs (Boone and Etchevers, 2001; Dutra et al., 2012; Wang et al., 2013; Cullather et al., 2014; Decharme et al., 2016; Born et al., 2019). These models have a limited number of layers and are based on simplified representations of the main processes affecting the SMB changes but usually do not have any explicit representation of snow metamorphism. However, they offer a good compromise between models of high complexity and simplified approaches or bulk-layer models for coupling with atmospheric models.

The snow module Explicit Snow (referred to hereafter as ES) implemented in the land surface model ORCHIDEE (Organising Carbon and Hydrology In Dynamic Ecosystems; Krinner et al., 2005; Chéruy et al., 2020) of the IPSL-CM ESM (Boucher et al., 2020) belongs to this third class of snow models. It has been successfully evaluated against observations in Col de Porte (French Alps) and in various sites in northern Eurasia (Wang et al., 2013). However, it has not yet been adapted to ice-covered areas. As a result, glaciers are considered bare soil in the current ORCHIDEE version, and over ice sheets, snow is handled with the atmospheric component of IPSL-CM in a very simplistic way. Recently, we made new developments to apply the ES model to glaciers and ice sheets, with a special focus on the GrIS. These developments meet two objectives. The first one is to treat snow-related processes in IPSL-CM in a more consistent way for all surface types. The second one is to compute the SMB, taking the main processes occurring within the snowpack into account. These developments also constitute a preliminary step for the subsequent use of the computed SMB as an interface between IPSL-CM and ice-sheet models. In the following, we will refer to ORCHIDEE-ICE when we mean the version of ORCHIDEE that includes these new developments and to ORCHIDEE for the former version of the model.

In this study, we evaluate the computation of SMB (and its components) in the ES model. As SMB is strongly dependent on the albedo, we also examine its sensitivity to various albedo parameterizations. To achieve this, we performed offline ORCHIDEE-ICE simulations and compared our results against model outputs from the polar-oriented regional atmospheric model MARv3.11.4 (Modèle Atmosphérique Régional; Fettweis et al., 2017) and the MODIS (MODerate resolution Imaging Spectroradiometer; Hall et al., 1995; Hall and Riggs, 2016) surface albedo retrievals. The paper is organized as follows. In Sect. 2, we provide an extensive description of the main characteristics of the original ES model as well as changes that have occurred since its early publication (Wang et al., 2013). The new developments made to apply ES to the GrIS are also presented in this section. Section 3 describes the experimental setup, and Sect. 4 provides a brief overview of the different datasets used for evaluation. The results are presented in Sects. 5 and 6 and discussed in Sect. 7.

2 Model description

2.1 Snow processes in the current ORCHIDEE-AR6 model

ORCHIDEE is the land surface component of the IPSL-CM Earth system model (Boucher et al., 2020; Chéruy et al., 2020) mainly developed at the French Institute Pierre Simon Laplace (IPSL). It computes both water and energy exchanges (SECHIBA module) between land surfaces and the atmosphere at a 30 min time step and includes carbon-related processes (STOMATE module). Within a given grid cell, land cover is represented as fractions of bare soils and vegetated areas described in terms of plant functional types (PFTs). Snow–vegetation interactions are not explicitly represented, and snow is evenly distributed among the various PFTs. Soil types are prescribed according to the USDA soil texture maps (Reynolds et al., 2000). The ORCHIDEE model can be run in offline mode, driven by atmospheric fields, or coupled with an atmospheric model. In the former ORCHIDEE version used for CMIP5 (Taylor et al., 2012), the snow scheme over glaciated surfaces was based on the bulk approach proposed by Chalita and Le Treut (1994). It consisted of a composite soil–snow model accounting for the thermal and radiative properties of snow cover (i.e. albedo and its variations with snow ageing). Snow was described as having a constant density (330 kg m⁻³), and melting occurred when the temperature exceeded 0 °C. Other processes such as water percolation and refreezing were ignored, although they directly impact the water budget. This means that all liquid water coming from melting snow left the snowpack as runoff.

For the CMIP6 exercise (Eyring et al., 2016), the bulk approach was replaced by the ES snow scheme, which was formerly adapted to the ORCHIDEE architecture (Wang et al., 2013) from a three-layer version of the ISBA-ES scheme (Interactions between Soil, Biosphere and Atmosphere-Explicit Snow scheme; Boone and Etchevers, 2001) developed at the French National Centre for Meteorological Research. The ES model is now used in the standard version of ORCHIDEE (version 2.0 onwards). However, it has not yet been considered for use over mountainous glaciers, which are treated as bare soils, nor over ice-sheet areas, which are currently handled by the LMDZ atmospheric model (Chéruy et al., 2020) with a very elementary snow scheme (i.e. a single-layer model, constant albedo, and thermal conductivity). In this section, we provide an extensive description of the snow model, including the main differences from the original ISBA-ES version (Wang et al., 2013). The new developments accounting for snow processes over ice-covered areas in the ORCHIDEE model are described in Sect. 2.2.

The ES model represents the snowpack as a one-dimensional physical system (vertical coordinate z). This means that all the lateral fluxes of mass and energy are ignored. The original version of this snowpack is discretized in three layers following the parameterization of Lynch-Stieglitz (1994), which sets the upper limits for the thickness of the first two layers to 5 and 50 cm respectively. This ensures the propagation of variations in the diurnal cycle of temperature and radiation and enables vertical heat and density gradients, which are assumed to be larger near the surface, to be resolved correctly. Each layer is described in terms of snow density, snow age, layer thickness, heat content, snow temperature, and liquid water content, with the first three variables being prognostic variables. Changes in snow mass are determined by the snowfall rate, snow melting, runoff at the base of the snowpack, and sublimation at the surface. In the absence of coupling with a dynamic ice-sheet model, snow mass at the surface of the ice sheet can be overestimated. Thus, to prevent excessive snow accumulation, we impose a maximum threshold of 3000 kg m⁻² beyond which snow is artificially removed. An overview of the organization of the different subroutines of the ES snowpack model is provided in Fig. 1. The description of the processes is given in the following subsections, and the list of model parameters is provided in Table A1 (Appendix A).

Figure 1 Flowchart of the new Explicit Snow scheme implemented in the ORCHIDEE-ICE model.

2.1.1 Surface processes

Energy balance

The evolution of the snowpack is primarily driven by the energy flux at the snow–atmosphere interface. A single energy balance is computed for all surface types coexisting in one grid cell. The surface energy flux (G_surf) available at the snow–atmosphere interface is computed from the energy balance equation:

$$\begin{array}{}\text{(1)}& G_{\mathrm{surf}}={\text{SW}}_{\mathrm{net}}+{\text{LW}}_{\mathrm{net}}-H_{\mathrm{L}}-H_{\mathrm{S}}+H_{\mathrm{rainfall}}.\end{array}$$

G_surf is computed positively when it warms the soil. SW_net and LW_net are the net shortwave and longwave radiation respectively, H_L is the latent heat flux, H_S is the sensible heat flux, and H_rainfall is the energy released by rainfall (see Eq. 14 in Boone and Etchevers, 2001). Equation (1) is used to compute the surface temperature (T_surf) of the grid cell at the next time step and provides the limit condition of the surface temperature at the snow–atmosphere interface for the calculation of the snow temperature profile.

Above snow-covered surfaces when T_surf is above the freezing temperature T₀ (273.15 K), the excess energy is first used to bring the snow temperature to T₀. A surface energy flux G_freezing associated with the freezing temperature T₀ can be computed using a similar formulation to Eq. (1). The difference between G_surf and G_freezing is converted to an additional temperature expressed as

$$\begin{array}{}\text{(2)}& T_{\mathrm{snow}}^{\mathrm{add}}=T_{\mathrm{surf}}-T_{\mathrm{0}}={\displaystyle \frac{G_{\mathrm{surf}}-G_{\mathrm{freezing}}}{C_{\mathrm{soil}}}}\mathrm{d}t.\end{array}$$

C_soil is the surface heat capacity of soil (J m⁻² K⁻¹) and is computed as the sum of heat capacities for snow-covered and snow-free surfaces (for both non-glaciated and glaciated areas), weighted by their respective grid cell fractions. For snow-covered surfaces, the specific heat capacity is defined as the product of snow density and the specific heat of ice (2106 J K⁻¹ kg⁻¹). If $T_{\mathrm{snow}}^{\mathrm{add}}$ is greater than (or equal to) the freezing temperature, the energy excess is used for melting snow, and G_surf is further set to G_freezing for energy conservation. If the new G_surf value is greater than the total heat content of the snowpack, snow is entirely melted and the excess energy is transferred to the underlying soil. The energy released by snowfall is accounted for in the snowpack scheme to update the snow heat content of the snowpack after a snowfall event.

Turbulent heat fluxes

The sensible (H_S) and latent heat (H_L) fluxes computed for each grid cell are given respectively by

$$\begin{array}{}\text{(3)}& {\displaystyle }{H}_{\mathrm{S}}={\mathit{\rho }}_{\mathrm{air}}{q}_{\mathrm{cdrag}}U\left(T_{\mathrm{surf}}-T_{\mathrm{air}}\right),\text{(4)}& {\displaystyle }{H}_{\mathrm{L}}=L_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{air}}{q}_{\mathrm{cdrag}}U\left(Q_{\mathrm{sat}}-Q_{\mathrm{air}}\right),\end{array}$$

where ρ_air is the air density, T_surf and T_atm are the surface and the 2 m atmospheric temperatures, Q_air and Q_sat are the air specific humidity at 2 m and the saturated specific humidity at the surface, L_s is the latent heat of sublimation (2.8345×10⁶ J kg⁻¹), U is the wind speed at 10 m, and q_cdrag is the drag coefficient computed as a function of the ice roughness length (z0_ice = 0.001 m) following the Monin–Obukhov turbulence theory (Monin and Obukhov, 1954) and the parameterizations of the eddy fluxes proposed by Louis (1979).

Snow sublimation

The amount of sublimation is simply deduced from the latent heat flux:

$$\begin{array}{}\text{(5)}& S_{\mathrm{snow}}={\displaystyle \frac{H_{\mathrm{L}}}{L_{\mathrm{s}}}}.\end{array}$$

Snow cover fraction

The snow cover fraction (F_snow) is derived from the formulation of Niu and Yang (2007), which has been shown to better represent the seasonal variation in the relationship between snow depth (Z_snow) and snow cover fraction thanks to its dependence on snow density:

$$\begin{array}{}\text{(6)}& F_{\mathrm{snow}}=\tanh \left({\displaystyle \frac{Z_{\mathrm{snow}}}{\mathrm{2.5}{z}_{\mathrm{0}\mathrm{g}}\times {\left(\frac{\langle {\mathit{\rho }}_{\mathrm{snow}}\rangle }{{\mathit{\rho }}_{\min }}\right)}^m}}\right),\end{array}$$

where 〈ρ_snow〉 is the snow density averaged over the total thickness of the snowpack; ρ_min is the minimum snow density (set to 50 kg m⁻³), that is, the density of fresh snow; z_0g is the ground roughness length (set to 0.01 m); and m (set to 1.0 in the present study) is an adjustable parameter.

Snow albedo

Compared to the early version presented in Wang et al. (2013), the albedo scheme has been modified, and snow albedo is now computed following the formulation of Chalita and Le Treut (1994):

$$\begin{array}{}\text{(7)}& {\mathit{\alpha }}_{\mathrm{snow}}=A_{\mathrm{aged}}+B_{\mathrm{dec}}\exp \left(-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{snow}}}{{\mathit{\tau }}_{\mathrm{dec}}}}\right),\end{array}$$

where A_aged represents the albedo of a snow-covered surface after snow ageing (old snow), and B_dec is defined so that the sum of A_aged and B_dec represents the albedo of fresh snow (i.e. maximum snow albedo). τ_dec is the time constant of the albedo decay and τ_snow is the snow age and is parameterized as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathit{\tau }}_{\mathrm{snow}}\left(t+\mathrm{d}t\right)& =\left[{\mathit{\tau }}_{\mathrm{snow}}\left(t\right)+\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{snow}}}{{\mathit{\tau }}_{\max }}}\right)\times \mathrm{d}t\right]\\ & \times \exp \left(-{\displaystyle \frac{P_{\mathrm{snow}}}{{\mathit{\delta }}_{\mathrm{c}}}}\right)+f_{\mathrm{age}},\end{array}\end{array}$$

where τ_max is the maximum snow age, P_snow is the amount of snowfall during the time interval dt, and δ_c is the critical value of solid precipitation necessary for reducing the snow age by a factor $\mathrm{1}/e$. As the ORCHIDEE time step is fixed to 30 min, the snow age is almost zero in a few time steps. In addition, the low surface air temperatures found in polar regions slow down the metamorphism. This effect is accounted for with the function f_age expressed as

$$\begin{array}{}\text{(9)}& {\displaystyle }{f}_{\mathrm{age}}=\left[{\displaystyle \frac{\left({\mathit{\tau }}_{\mathrm{snow}}\left(t\right)+\left(\mathrm{1}-\frac{{\mathit{\tau }}_{\mathrm{snow}}}{{\mathit{\tau }}_{\max }}\right)\times \mathrm{d}t\right)\times \exp \left(-\frac{P_{\mathrm{snow}}}{{\mathit{\delta }}_{\mathrm{c}}}\right)-{\mathit{\tau }}_{\mathrm{snow}}\left(t\right)}{\mathrm{1}+g_{\mathrm{temp}}\left(T_{\mathrm{surf}}\right)}}\right]\text{(10)}& {\displaystyle }{g}_{\mathrm{temp}}\left(T_{\mathrm{surf}}\right)={\left[{\displaystyle \frac{max(T_{\mathrm{0}}-T_{\mathrm{surf}},\mathrm{0})}{{\mathit{\omega }}_{\mathrm{1}}}}\right]}^{{\mathit{\omega }}_{\mathrm{2}}},\end{array}$$

where ω₁ and ω₂ are tuning constants. The albedo is computed for the visible and near-infrared spectral bands. However, to compute the upward shortwave radiation, an arithmetic mean between the visible and the near-infrared albedo is considered.

A single energy balance is computed for all surface types, but the albedo is weighted by the different fractions of PFTs and glaciated areas and by the snow-covered and snow-free fractions. As a result, the surface albedo (α) of the grid cell is computed as the sum of snow-free albedo (α_snow-free) and snow-covered albedo (α_snow) weighted by the fractional area of the grid cell covered by snow F_snow (snow-covered fraction hereafter):

$$\begin{array}{}\text{(11a)}& \mathit{\alpha }=F_{\mathrm{snow}}\times {\mathit{\alpha }}_{\mathrm{snow}}+\left(\mathrm{1}-F_{\mathrm{snow}}\right)\times {\mathit{\alpha }}_{\text{snow-free}},\end{array}$$

with

$$\begin{array}{}\text{(11b)}& {\mathit{\alpha }}_{\mathrm{snow}}=f_{\mathrm{ice}}\times {\mathit{\alpha }}_{\mathrm{snow}}^{\mathrm{ice}}+\sum _{\mathrm{PFT}}{f}_{\mathrm{PFT},i}\times {\mathit{\alpha }}_{\mathrm{snow}}^{\mathrm{PFT},i}\end{array}$$

and

$$\begin{array}{}\text{(11c)}& {\mathit{\alpha }}_{\text{snow-free}}=f_{\mathrm{ice}}\times {\mathit{\alpha }}_{\text{snow-free}}^{\mathrm{ice}}+\sum _{\mathrm{PFT}}{f}_{\mathrm{PFT},i}\times {\mathit{\alpha }}_{\text{snow-free}}^{\mathrm{PFT},i}.\end{array}$$

f_ice and f_PFT,i are the grid cell fractions of ice-covered areas and the ith PFT respectively; ${\mathit{\alpha }}_{\mathrm{snow}}^{\mathrm{ice}}$ (${\mathit{\alpha }}_{\text{snow-free}}^{\mathrm{ice}}$) and ${\mathit{\alpha }}_{\mathrm{snow}}^{\mathrm{PFT},i}$ (${\mathit{\alpha }}_{\text{snow-free}}^{\mathrm{PFT},i}$) are the corresponding snow albedo (snow-free albedo) values.

Over the GrIS, α_snow-free is given by the albedo of bare ice, prescribed to 0.6 and 0.2 for the visible and near-infrared wavelengths respectively. At the margins of the GrIS, some grid points may be only partially covered by snow or ice or even become totally snow-free during the melting season. It is therefore important to take these different features into account to correctly compute the surface albedo of the GrIS.

2.1.2 Internal processes

When snow falls on a snow-free surface, a new snowpack is generated provided that the ground temperature is below or equal to the freezing point. The snow mass and the heat content of the snowfall are initially distributed evenly within the three layers. The snow density is the same for the three layers and is given by the density of the snowfall computed as a function of wind speed and surface air temperature (Pahaut, 1976). When snowfall occurs over an existing snowpack, fresh snow is added to the upper layer providing that the snowfall thickness is greater than the critical threshold δ_c (see Eq. 8). The snow thickness, density, and heat content are then modified in this layer. However, as the number of snow layers is kept fixed, redistribution of mass and heat content within the layers is required when snow depth changes, but the total snow mass and heat content are conserved.

Heat conduction

Solar absorption is not accounted for in the snow model. All incoming solar energy is therefore deposited at the snow surface and distributed to deeper layers through heat conduction. The heat conduction from the surface to the bottom of the snowpack is described by a vertical-diffusion equation relating the temporal evolution of the snow temperature in the snowpack at depth z and the divergence of the snow heat flux F_C and is solved using an implicit numerical scheme.

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial T_{\mathrm{snow}}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{C_{\mathrm{snow}}}}\cdot {\displaystyle \frac{\partial F_C}{\partial z}}\text{(13)}& {\displaystyle }{F}_C=-{\mathrm{\Lambda }}_{\mathrm{s}}{\displaystyle \frac{\partial T_{\mathrm{snow}}}{\partial z}},\end{array}$$

with C_snow (J m⁻² K⁻¹), Λ_s, and T_snow being the snow heat capacity, the snow thermal conductivity (W m⁻¹ K⁻¹), and the snow temperature respectively.

At the snow–atmosphere interface, the boundary condition is given by the energy balance equation (F_C=G_surf) and is used in the ORCHIDEE model to compute the surface temperature.

Along with the thermal gradient, a water vapour diffusive flux takes place from the warmer to the colder parts of the snowpack, and sublimation or condensation may occur in the pore spaces depending on the water vapour saturation pressure. This process is particularly significant in the Arctic because of strong temperature gradients between soils and atmosphere, and it is largely responsible for snow metamorphism. While it is explicitly accounted for in detailed snow models, in Explicit Snow, the effect of water vapour diffusion and phase changes is parameterized through the thermal conductivity (Sun et al., 1999). The effective thermal conductivity (Λ_eff) is thus expressed as the sum of empirical formulations for snow thermal conductivity (Λ_cond) and thermal conductivity from vapour transport (Λ_vap), with

$$\begin{array}{}\text{(14)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{cond}}^i=a_{\mathit{\lambda }}+b_{\mathit{\lambda }}{\mathit{\rho }}_{\mathrm{snow}}^{i^{\mathrm{2}}}\text{(15)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{vap}}^i=\left(a_{\mathit{\lambda }v}+{\displaystyle \frac{b_{\mathit{\lambda }v}}{c_{\mathit{\lambda }v}+T_{\mathrm{snow}}^i}}\right){\displaystyle \frac{P_{\mathrm{0}}}{P}},\end{array}$$

where a_λ=0.02 W m⁻¹ K⁻¹, $b_{\mathit{\lambda }}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{6}}$ W m⁵ K⁻¹ kg⁻² (Anderson, 1976), $a_{\mathit{\lambda }v}=-\mathrm{0.06023}$ W m⁻¹ K⁻¹, $b_{\mathit{\lambda }v}=-\mathrm{2.5425}$ W m⁻¹, and $c_{\mathit{\lambda }v}=-\mathrm{289.99}$ K (Yen, 1981). P is the atmospheric pressure (in hPa) and P₀=1000 hPa. The superscripts i denote the ith layer.

Heat content

The heat content is computed using the following equation:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{H}_{\mathrm{snow}}^i& =D_{\mathrm{snow}}^i\left[C_{\mathrm{snow}}^{v,i}\left(T_{\mathrm{snow}}^i-T_{\mathrm{f}}\right)-L_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{snow}}^i\right]\\ & +L_{\mathrm{f}}{\mathit{\rho }}_{\mathrm{water}}{W}_{\mathrm{liq}}^i,\end{array}\end{array}$$

where L_f is the latent heat of fusion, and ρ_water is the water density. $H_{\mathrm{snow}}^i$, $W_{\mathrm{liq}}^i$, $D_{\mathrm{snow}}^i$, ${\mathit{\rho }}_{\mathrm{snow}}^i$, and $C_{\mathrm{snow}}^{v,i}$ are the heat and liquid contents, the depth, the density, and the mean volumetric heat capacity (J K⁻¹ m⁻³) of the ith layer.

After heat redistribution within the snowpack, snow temperature is diagnosed using Eq. (16), assuming no liquid water in the snowpack. If the snow temperature exceeds the freezing point, the liquid content in each layer is then diagnosed from the snow temperature and heat content of the layer, and the temperature is then reset to the freezing point.

Compaction

The total snow depth decreases as density increases. Changes in density occur as a result of the weight of the overlying snow layers and under the influence of snow metamorphism. The local rate of density change in the ith layer is derived from Anderson (1976):

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{snow}}^i}}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{snow}}^i}{\partial t}}& ={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{snow}}^i}{{\mathit{\eta }}_{\mathrm{snow}}^i\left(T_{\mathrm{snow}}^i,{\mathit{\rho }}_{\mathrm{snow}}^i\right)}}\\ & +{\mathit{\psi }}_{\mathrm{snow}}^i\left(T_{\mathrm{snow}}^i,{\mathit{\rho }}_{\mathrm{snow}}^i\right).\end{array}\end{array}$$

The first term on the right-hand side represents the compaction due to snow load, with ${\mathit{\sigma }}_{\mathrm{snow}}^i$ (Pa) being the pressure of the overlying snow and ${\mathit{\eta }}_{\mathrm{snow}}^i$ the snow viscosity.

$${\mathit{\sigma }}_{\mathrm{snow}}^i=g\times M_{\mathrm{snow}}^i,$$

where g is the gravitational constant (m s⁻²) and $M_{\mathrm{snow}}^i$ the cumulative snow mass (kg m⁻²).

The viscosity (in Pa s) is expressed as a function of snow temperature and density (Mellor, 1964; Kojima, 1967):

$$\begin{array}{}\text{(18)}& {\mathit{\eta }}_{\mathrm{snow}}^i={\mathit{\eta }}_{\mathrm{0}}\exp \left[a_{\mathit{\eta }}\left(T_{\mathrm{f}}-T_{\mathrm{snow}}^i\right)+b_{\mathit{\eta }}{\mathit{\rho }}_{\mathrm{snow}}^i\right],\end{array}$$

with ${\mathit{\eta }}_{\mathrm{0}}=\mathrm{3.7}\times {\mathrm{10}}^{\mathrm{7}}$ Pa s, $a_{\mathit{\eta }}=\mathrm{8.1}\times {\mathrm{10}}^{-\mathrm{2}}$ K⁻¹, and $b_{\mathit{\eta }}=\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{2}}$ m³ kg⁻¹.

The second term on the right-hand side of Eq. (17) parameterizes the effect of metamorphism, which is significant for newly fallen snow.

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\mathit{\psi }}_{\mathrm{snow}}^i& =a_{\mathit{\psi }}\exp \left[-b_{\mathit{\psi }}\cdot \left(T_{\mathrm{f}}-T_{\mathrm{snow}}^i\right)\right.\\ & \left.-c_{\mathit{\psi }}\cdot max\left(\mathrm{0},{\mathit{\rho }}_{\mathrm{snow}}^i-{\mathit{\rho }}_{\mathit{\psi }}\right)\right]\end{array}\end{array}$$

The values of the parameters are the following: $a_{\mathit{\psi }}=\mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻¹, $b_{\mathit{\psi }}=\mathrm{4.2}\times {\mathrm{10}}^{-\mathrm{2}}$ K⁻¹, c_ψ=460 m³ kg⁻¹, and ρ_ψ=150 kg m⁻³.

In the model, density changes due to compaction are allowed as long as density remains below a threshold, fixed at 750 kg m⁻³. This value was chosen because compaction becomes slower above densities between 550 and 800 kg m⁻³ due to the progressive disappearance of air spaces between the snow particles (Maeno and Ebinuma, 1983). A critical value of 730 kg m⁻³ has even been proposed by Maeno (1978). Compaction does not affect the total mass and the heat content of the snowpack but changes the layer thicknesses. The distribution of snow heat within the layers must therefore be updated using Eq. (16).

Vertical temperature profile

The snow temperature profile resulting from heat redistribution is then computed by solving the heat diffusion equation using an implicit numerical scheme similar to that used for heat diffusion in the soil. The vertical temperature profile within the snowpack is expressed as

$$\begin{array}{}\text{(20)}& T_{\mathrm{snow}}^{\mathrm{1}}=\left[{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{snow}}\cdot C_{\mathrm{gr}\mathrm{\_}\mathrm{snow}}+\left(T_{\mathrm{surf}}+T_{\mathrm{snow}}^{\mathrm{add}}\right)}{\mathrm{1}+{\mathit{\lambda }}_{\mathrm{snow}}\left(\mathrm{1}-D_{\mathrm{gr}\mathrm{\_}\mathrm{snow}}\right)}}\right]\end{array}$$

for the first layer and as

$$\begin{array}{}\text{(21)}& T_{\mathrm{snow}}^{i+\mathrm{1}}=C_{\mathrm{gr}\mathrm{\_}\mathrm{snow}}+D_{\mathrm{gr}\mathrm{\_}\mathrm{snow}}\cdot T_{\mathrm{snow}}^i\end{array}$$

for the deeper layers (i>1), where λ_snow, C_{gr_snow} and D_{gr_snow} are coefficients resulting from the resolution of the numerical scheme and depend on the snow heat capacity, thermal conductivity, and characteristics of the vertical discretization. The numerical scheme is similar to the one presented in Wang et al. (2016; see appendix A therein) in which the temperature at the interface between two layers is calculated as a linear interpolation according to the two nearest nodes (middle of the layers). Diffusion therefore takes place downward and upward.

Melting and refreezing processes

If meltwater is produced at the surface, it may remain in a liquid state in the uppermost layer or penetrate the next layer, where it can remain or refreeze as long as the maximum water-holding capacity is not reached; otherwise, it penetrates the lower layers.

The evolution of liquid water in each layer is controlled by the energy available to induce phase changes and by the maximum water-holding capacity. In the ith layer, the energy used for melting snow ($E_{\mathrm{snow}}^i$) is expressed as

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{E}_{\mathrm{snow}}^i& =min\left(C_{\mathrm{snow}}^{v,i}{D}_{\mathrm{snow}}^i\times max\left(\mathrm{0},T_{\mathrm{snow}}^i-T_{\mathrm{f}}\right)\right.,\\ & \left.max\left(\mathrm{0},D_{\mathrm{swe}}^i-W_{\mathrm{liq}}^i\right)\times L_{\mathrm{f}}{\mathit{\rho }}_{\mathrm{water}}\right),\end{array}\end{array}$$

where $D_{\mathrm{swe}}^i$ is the snow water equivalent in the ith layer. The first term represents the available energy for phase change in the ith layer, and the second term corresponds to the energy required to entirely melt the snow mass that has not been transformed into liquid water. The maximum water-holding capacity is taken from Anderson (1976):

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{W}_{\max }^i& =\left[r_{\min }+\left(r_{\max }-r_{\min }\right)\cdot max\left(\mathrm{0},{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{t}}-{\mathit{\rho }}_{\mathrm{snow}}^i}{{\mathit{\rho }}_{\mathrm{t}}}}\right)\right]\\ & \cdot {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{snow}}^i}{{\mathit{\rho }}_{\mathrm{w}}}}\cdot D_{\mathrm{snow}}^i,\end{array}\end{array}$$

with r_min=0.03, r_max=0.10, and ρ_t=200 kg m⁻³.

Runoff (S_melt) is computed as the sum of meltwater produced at the surface and the total liquid water that has percolated down to the bottom layer and that exceeds $W_{\max }^{\mathrm{bottom}}$. It is thus simply given by

$$\begin{array}{}\text{(24)}& M_{\mathrm{snow}}={\displaystyle \frac{\sum _i{E}_{\mathrm{snow}}^i}{L_{\mathrm{f}}}}.\end{array}$$

At each time step, changes in layer thickness, density, and liquid water content in each layer are updated, as well as changes in snow temperature due to melting or refreezing. In case of complete snow melting, the energy excess that has not been used for phase changes is used to warm the underlying ground.

2.2 New developments

2.2.1 New snow-layering scheme

As mentioned in Sect. 1, snow models of intermediate complexity are a good compromise between detailed snow models and single-layer models. They are designed to be implemented in ESMs and as such should not require excessive computational time. Although their vertical resolution is generally limited to five layers at most (Cristea et al., 2022), several studies reported that snow models of intermediate complexity considerably improve the representation of basic features of the snowpack and reduce biases in surface temperature when they are compared to single-layer models (Lynch-Stieglitz, 1994; Boone and Etchevers, 2001; Dutra et al., 2012; Wang et al., 2013). Despite this good performance, increasing the number of snow layers (with finer layers near the surface or near the snow–ice interface) is expected to improve the modelled heat conduction within the snowpack, the simulated temperature at the snow–ice interface, and subsequently the vertical temperature profile in the ice and eventually the simulated SMB (Cristea et al., 2022). We therefore increased the number of snow layers from 3 to 12, following the layering scheme proposed by Decharme et al. (2016) for ISBA-ES, in which the new layering scheme is defined as

$$\begin{array}{}\text{(25)}& \left\{\begin{array}{ll}{D}_{\mathrm{snow}}^i& =min\left({\mathit{\delta }}_i,\frac{Z_{\mathrm{snow}}}{\mathrm{12}}\right)\text{ for }i\le \mathrm{5}\text{ or }i\ge \mathrm{9}\\ & \\ D_{\mathrm{snow}}^{\mathrm{6}}& =\mathrm{0.3}{d}_{\mathrm{r}}-min\left(\mathrm{0},\mathrm{0.3}{d}_{\mathrm{r}}-D_{\mathrm{snow}}^{\mathrm{5}}\right)\\ & \\ D_{\mathrm{snow}}^{\mathrm{7}}& =\mathrm{0.4}{d}_{\mathrm{r}}+min\left(\mathrm{0},\mathrm{0.3}{d}_{\mathrm{r}}-D_{\mathrm{snow}}^{\mathrm{5}}\right)\\ & -min\left(\mathrm{0},\mathrm{0.3}{d}_{\mathrm{r}}-D_{\mathrm{snow}}^{\mathrm{9}}\right)\\ & \\ D_{\mathrm{snow}}^{\mathrm{8}}& =\mathrm{0.3}{d}_{\mathrm{r}}-min\left(\mathrm{0},\mathrm{0.3}{d}_{\mathrm{r}}-D_{\mathrm{snow}}^{\mathrm{9}}\right)\\ & \\ d_{\mathrm{r}}& =Z_{\mathrm{snow}}-\sum _{i=\mathrm{1}}^{\mathrm{5}}{D}_{\mathrm{snow}}^i-\sum _{i=\mathrm{9}}^{\mathrm{12}}{D}_{\mathrm{snow}}^i\end{array}.\right.\end{array}$$

The δ_i values correspond to the maximum widths of the layers 1 to 5 and 9 to 12 and are fixed to δ₁=0.01 m, δ₂=0.05 m, δ₃=0.15 m, ${\mathit{\delta }}_{\mathrm{4}}={\mathit{\delta }}_{\mathrm{10}}=\mathrm{0.5}$ m, ${\mathit{\delta }}_{\mathrm{5}}={\mathit{\delta }}_{\mathrm{9}}=\mathrm{1}$ m, δ₁₁=0.1 m, and δ₁₂=0.02 m. For very thin snowpacks ($Z_{\mathrm{snow}}\le Z_{\mathrm{thin}}=\mathrm{0.1}$ m), each layer has the same thickness $\frac{Z_{\mathrm{thin}}}{\mathrm{12}}$ The layer thicknesses are updated at each time step if the first two layers (i=1, 2) or the bottom layer (i=12) become too thin $\left(\text{less than }{D}_{\mathrm{snow}}^i=\mathrm{0.5}\times min\left({\mathit{\delta }}_i,\frac{Z_{\mathrm{snow}}}{\mathrm{12}}\right)\right)$ or too thick $\left(\text{larger than }{D}_{\mathrm{snow}}^i=\mathrm{1.5}\times min\left({\mathit{\delta }}_i,\frac{Z_{\mathrm{snow}}}{\mathrm{12}}\right)\right)$. In that case, the snow mass and heat content are redistributed according to the new layering scheme. Otherwise, the layer thicknesses at the current time step are kept to their previous values (i.e. at the previous time step). This allows us to maintain the density and thermal conductivity of fresh snow as long as the depth has not changed too much. This enables the model to work more closely with more complex models in which new snow layers are associated with a new snowfall event.

2.2.2 Implementation of ice layers

In the case when the snow mass melts completely, ice melting occurs if the available energy is sufficient and contributes to runoff. To account for the presence of ice below the snow layers, we implemented a new module in ORCHIDEE to compute the heat diffusion and the vertical temperature distribution in the ice, as well as the potential ice melting. This module works in a similar way as the ES model and only accounts for vertical fluxes. The ice reservoir is discretized into eight layers whose maximum thicknesses are fixed at 0.01, 0.05, 0.15, 0.5, 1, 5, 10, and 50 m. A finer vertical spacing is imposed for the upper layers to better resolve heat conduction at the snow–ice or atmosphere–ice interface. The large thickness of the bottom layer allows it to have an almost-constant temperature throughout the year, as has been observed at a few tens of metres in depth (Patterson, 1994). Ice layers are only implemented above an icy soil type. If the icy soil is predominant in a given grid cell, then the entire surface corresponding to this grid point will be considered icy.

In the absence of a dynamic ice model that transports ice from the interior of the ice sheet (or glacier) to the edges, the total ice mass may disappear entirely in the ablation zones, especially in long-term simulations. To avoid such situations, ice is considered an infinite reservoir: melting ice contributes to runoff but, at each time step, the amount of ice melted in the upper layers is counterbalanced by ice added at the base, and the layer thicknesses are kept fixed to their initial value.

The vertical distribution of temperature is determined using the same numerical scheme as that for the snowpack. If the snow is still present over the icy soil, the temperature in the top ice layer is given by the temperature of the bottom snow layer computed using Eq. (21). If the snow has completely melted, the temperature in the first ice layer is given by an expression similar to Eq. (20):

$$\begin{array}{}\text{(26)}& T_{\mathrm{ice}}^{\mathrm{1}}=\left[{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{ice}}\cdot C_{\mathrm{gr}\mathrm{\_}\mathrm{ice}}+\left(T_{\mathrm{surf}}+T_{\mathrm{snow}}^{\mathrm{add}}\right)}{\mathrm{1}+{\mathit{\lambda }}_{\mathrm{ice}}\left(\mathrm{1}-D_{\mathrm{gr}\mathrm{\_}\mathrm{ice}}\right)}}\right].\end{array}$$

For the deeper layers, the ice temperature is expressed as follows:

$$\begin{array}{}\text{(27)}& T_{\mathrm{ice}}^{i+\mathrm{1}}=C_{\mathrm{gr}\mathrm{\_}\mathrm{ice}}+D_{\mathrm{gr}\mathrm{\_}\mathrm{ice}}\cdot T_{\mathrm{ice}}^i.\end{array}$$

Similarly to the snow coefficients (see Eqs. 20 and 21), λ_ice, C_{gr_ice}, and D_{gr_ice} depend on the vertical discretization and the thermal properties of the ice. The formulations of the heat capacity (C_ice) and thermal conductivity (Λ_ice) of the ice have been taken from those used in the GRISLI ice-sheet model (Yen, 1981) and are given by

$$\begin{array}{}\text{(28)}& {\displaystyle }{C}_{\mathrm{ice}}={\mathit{\rho }}_{\mathrm{ice}}\left(a_{\mathrm{ci}}+b_{\mathrm{ci}}\left(T_{\mathrm{ice}}-T_{\mathrm{0}}\right)\right)\text{(29)}& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{ice}}=a_{\mathit{\lambda }i}\exp \left(b_{\mathit{\lambda }i}\times T_{\mathrm{0}}\right),\end{array}$$

where T_ice is the ice temperature, a_ci=2115.3 J K⁻¹ kg⁻¹, b_ci=7.79293 J K⁻² kg⁻¹, a_λi=6.727 W m⁻¹ K⁻¹, and $b_{\mathit{\lambda }i}=-\mathrm{0.041}$ K⁻¹.

A major difference between the hydrology of snow and ice layers lies in the fact that ice is considered an impermeable medium. Hence, liquid water coming from melting ice is treated as running off instantaneously with no possibility of refreezing. As a result, when the ice temperature is above the melting point, the energy available for phase change in the ith ice layer (J m⁻²) is given by

$$\begin{array}{}\text{(30)}& E_{\mathrm{ice}}^i=C_{\mathrm{ice}}^i\left(T_{\mathrm{ice}}^i-T_{\mathrm{0}}\right)D_{\mathrm{ice}}^i.\end{array}$$

Similarly to S_melt (Eq. 24), the total amount of ice melt is given by

$$\begin{array}{}\text{(31)}& M_{\mathrm{ice}}={\displaystyle \frac{\sum _i{E}_{\mathrm{ice}}^i}{L_{\mathrm{f}}}},\end{array}$$

and the runoff is computed as the sum of M_snow and M_ice. Given the fact that snow drift is ignored, the surface mass balance is computed as

$$\begin{array}{}\text{(32)}& \text{SMB}=P_{\mathrm{snow}}+P_{\mathrm{rain}}-M_{\mathrm{snow}}-M_{\mathrm{ice}}-S_{\mathrm{snow}}.\end{array}$$

2.2.3 Other processes in the new ES model

Another modification made to the ES module concerns the inclusion of rainwater percolation within the snowpack, which may refreeze at depth as long as the maximum water-holding capacity is not exceeded. In case of rain-on-snow events, we also enhanced the snow ageing by a factor of 2 ($f_{\mathrm{age}}=f_{\mathrm{age}}\times \mathrm{2}$). Although it sounds somewhat arbitrary, we introduced this parameterization into the model to account for the effect of such events on metamorphism and densification (Marshall et al., 1999), thereby lowering the albedo (Yang et al., 2023).

The snow thermal conductivity has been modified to follow a similar formulation to that used in the ISBA-ES model (Decharme et al., 2016) and the CROCUS model (Vionnet et al., 2012) and proposed earlier by Yen (1981). Therefore, the effective thermal conductivity in the ith layer now reads

$$\begin{array}{}\text{(33)}& {\mathrm{\Lambda }}_{\mathrm{eff}}^i=\left(a_{\mathit{\lambda }v}+{\displaystyle \frac{b_{\mathit{\lambda }v}}{c_{\mathit{\lambda }v}+T_{\mathrm{snow}}^i}}\right){\displaystyle \frac{P_{\mathrm{0}}}{P}}+{\mathrm{\Lambda }}_{\mathrm{ice}}{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}^i}{{\mathit{\rho }}_{\mathrm{w}}}}\right)}^{\mathrm{1.88}}.\end{array}$$

The first term on the right-hand side that parameterizes the water vapour diffusion effects (${\mathrm{\Lambda }}_{\mathrm{vap}}^i$) remains unchanged (see Eq. 15). The second term replaces Eq. (14) used in the previous ES version (Wang et al., 2013) and corresponds to the new formulation of the snow thermal conductivity (${\mathrm{\Lambda }}_{\mathrm{cond}}^i$). Here, the ice thermal conductivity (Λ_ice) differs from the value found in Decharme et al. (2016) and is given by Eq. (29).

Besides the new snow-layering scheme and the changes mentioned in this section, all the other processes simulated in the new ES module are treated in the same way as in the three-layer version.

3 Experimental setup

3.1 Forcing by the regional atmospheric model MAR

The ORCHIDEE-ICE simulations presented in this paper were driven by the atmospheric outputs of the regional atmospheric model MAR (Fettweis et al., 2017). This approach was motivated by the fact that MAR was initially developed for polar regions (Gallée and Schayes, 1994). Moreover, it is coupled to a land surface scheme, SISVAT (Soil Ice Snow Vegetation Atmosphere Transfer; De Ridder and Schayes, 1997), that includes a physically based snowpack model derived from the multi-layered snow model CROCUS (Brun et al., 1989, 1992). As such, MAR has been extensively used to simulate the present-day climate and surface mass balance of the GrIS and compares well to reanalyses and available data of SMB measurements (e.g. Fettweis et al., 2017, 2020; Franco et al., 2012; Montgomery et al., 2020; Delhasse et al., 2020). Therefore, the use of atmospheric forcings from MAR offers a good opportunity to assess the performance of our snow model at simulating the SMB and ablation-related processes.

The MAR simulations (1960–2019) were run at a 20 km × 20 km resolution. Here, we use version 3.11.4, identical to version 3.11.5 for the Greenland ice sheet (Smith et al., 2023). MAR was forced every 6 h at its lateral boundaries by the meteorological fields (temperature, humidity, wind, and pressure) coming from the ERA-40 (1960–1978, Uppala et al., 2005) and the ERA-Interim (1979–2019, Dee et al., 2011) reanalyses from the European Centre for Medium-Range Weather Forecasts (ECMWF). Sea surface temperatures and sea ice cover, also coming from ECMWF reanalyses, were 6 h prescribed.

3.2 The ORCHIDEE-ICE simulations

The ORCHIDEE-ICE simulations are run at a 30 min time step with the same spatial resolution as the MAR outputs (20 km × 20 km). The integration domain covers the whole of Greenland. ORCHIDEE-ICE is forced every 3 h by the downward shortwave and longwave radiation, the surface air temperatures, and the specific humidity (all at 2 m); the wind speed (at 10 m); the surface pressure; and the precipitation rate (split between rainfall and snowfall). Simulations are performed over the 1995–2019 period. The first 5 years (1995 to 1999) are used for the initialization of the snowpack and are not included in the analysis. However, to obtain reasonable thermal conditions within the ice layers, longer time integration is required. Thus, we performed a preliminary spin-up experiment over the same 25 years to infer an initial vertical temperature profile for the subsequent ORCHIDEE-ICE simulations.

The name and the characteristics of the different experiments presented in this paper are summarized in Table 1. Using the experimental design described above, we first ran the ES model with 3 and 12 snow layers (STD-3L and STD-12L experiments respectively) to evaluate the added value of the new layering scheme. These experiments were carried out with the albedo parameters used in the CMIP6 ORCHIDEE version (Chéruy et al., 2020) and referred to hereafter as the standard snow albedo parameters.

Due to the strong sensitivity of the SMB to the albedo, we also conducted two additional experiments with modified values of the albedo parameters. In the ASIM-12L experiment, we used the parameters inferred from the approach of Raoult et al. (2023). This latter was based on a data assimilation experiment using the MODIS retrievals. The main goal of their study was to optimize the albedo parameters so as to improve the albedo for the ice sheet as a whole, while giving an extra weight to the edges where the greatest amount of runoff is produced. In doing this, they also succeeded in improving the model–data fit between the ORCHIDEE albedo and MODIS retrievals over the whole GrIS and in reducing the root-mean-square error (RMSE) by ∼25 %. However, their work was done with a previous version of the ORCHIDEE-ICE model with only three snow layers and in which the ice layers were not implemented. Instead, ice was mimicked by a soil type whose porosity and volumetric water content were set to 98 % to simulate soil filled with frozen water.

The logical follow-up to the work of Raoult et al. (2023) would have been to apply the optimization algorithm to the new version of ORCHIDEE-ICE. Since this approach is highly time-consuming, it has not yet been carried out, but it will be the focus of further investigations. Therefore, using the new ORCHIDEE-ICE model version, we adopted a manual tuning approach (i.e. trial-and-error method) to adjust the albedo parameters (OPT-12L experiment). This procedure consists in (1) changing the parameter values, the new value being taken from the range reported in Table 1; (2) running the model with the new parameter values; (3) evaluating the model performance (in terms of SMB and its components) using statistical criteria (e.g. RMSE between MAR and ORCHIDEE-ICE); and (4) repeating steps (1) to (3) until an acceptable calibration is obtained (i.e. acceptable values of SMB, runoff, refreezing, and sublimation).

Finally, to assess the impact of the climatic fields used as inputs of ORCHIDEE-ICE, we performed another experiment (ERA5-12L experiment) by forcing the model with the ERA5 reanalysis (Hersbach et al., 2020) and using the same albedo parameters as in OPT-12L experiment.

Table 1 List of the ORCHIDEE-ICE experiments (first column) with values chosen for the different albedo parameters (standard albedo parameters for STD-3L and STD-12L; optimized albedo parameters inferred from Raoult et al., 2023, for ASIM-12L; and manually tuned parameters for OPT-12L and ERA-12L). Values in brackets indicate the range of values considered in the manual tuning approach for each parameter. Bold formatting indicates the parameter values used in the OPT-12L (manual tuning) experiment.

4 Methodology for the model performance evaluation

4.1 Comparison with MAR outputs

Our first objective is to assess the performance of the ORCHIDEE ICE model in representing the GrIS SMB. The period under study spans the 2000–2019 period. As mentioned in Sect. 3, MAR has revealed good ability at simulating the SMB of present-day Greenland when compared to observational data. Therefore, at the scale of the entire GrIS, our evaluation is made with respect to the MAR outputs (Figs. 2a–5a). In all simulations presented in this paper except ERA5-12L, the forcing fields of the ORCHIDEE-ICE model are provided by MAR outputs. These include solid and liquid precipitation, which constitute the accumulation (and the climatic) component of the SMB. Using the MAR forcing, our analysis of the ability of ORCHIDEE-ICE to reproduce ablation processes (runoff and sublimation) is made easier and is not biased by the use of another forcing.

4.2 Comparison with available data

In this study, we compared the albedo computed in ORCHIDEE-ICE with satellite-derived estimates of daily albedo. We used collection 6 from the MOD10A1 product (Hall et al., 1995) retrieved from the NASA spaceborne sensor MODIS. We chose this product because it has good spatiotemporal coverage over snow-covered areas. It is also one of the best-performing products in terms of comparison with in situ data (Urraca et al., 2022, 2023). Moreover, while studies based on the previous collection, Collection 5, reported deficiencies at latitudes higher than 70° N (Alexander et al., 2014), substantial improvements have been made to Collection 6 using all available observations for the acquisition period versus only four observations per day in Collection 5 (https://lpdaac.usgs.gov/products/mcd43d11v006/, last access: 22 January 2024). As a result, better-quality retrievals are obtained at high latitudes despite a slight negative bias (Urraca et al., 2022). To avoid inaccuracies in retrieved data due to the presence of clouds or aircraft condensation trails, the MOD10A1 albedo product used in this study was further processed by Box et al. (2017): data have been de-noised, gap-filled, corrected for the sun-angle bias, and validated using daily ground albedo values from the PROMICE (Programme for Monitoring of the Greenland ice sheet, Fausto et al., 2021) and GC-net automatic weather stations (Box et al., 2017).

We aggregated the albedo data (500 m × 500 m) onto the MAR grid to make the comparison between MODIS data and the ORCHIDEE-ICE outputs. In this study, we used the albedo data covering the 2000–2017 period because data for the years 2018 and 2019 were undefined. The resulting dataset may be used to calibrate the mean ORCHIDEE-ICE albedo, computed as the mean between the visible (from 0.4 to 0.7 µm) and near-infrared (from 0.7 to 2.5 µm) bands (see Sect. 2).

As in Fettweis et al. (2020), we also evaluated the modelled SMB using the Machguth et al. (2016) SMB database. Daily outputs are used here over 2000–2019. Modelled SMBs were linearly interpolated to the measurement point location and corrected for the elevation difference between the MAR native topography at 20 km and the one provided in the SMB database. This was done using space-varying SMB–elevation gradients as proposed by Franco et al. (2012) and Noël et al. (2016). Finally, measurements not included in the 2000–2019 period and records located outside the 20 km MAR ice mask are discarded from the evaluation.

4.3 Statistical metrics

To evaluate our model performance, we used statistical metrics.

The root-mean-square error (RMSE) has been computed using the monthly mean variables averaged over 2000–2019 for the SMB and its components and over 2000–2017 for the albedo. It is defined as

$$\begin{array}{}\text{(34)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N{\left(X_{\mathrm{OR}}\left(i\right)-X_{\mathrm{MAR}}\left(i\right)\right)}^{\mathrm{2}}},\end{array}$$

where X_OR(i) and X_MAR(i) represent the ORCHIDEE-ICE and MAR variables respectively at each grid point i, N is the number of unmasked grid points (i.e. related only to the ice-covered area), and i stands for the ith grid point. The RMSE is a metric widely used to compare different models, but it has some shortcomings related to the fact that higher weights are given to larger errors. We there used additional statistical criteria to provide a more in-depth picture of our analysis. We computed the spatial RMSE (SRMSE), which gives a measure of the quadratic difference averaged over time between values simulated by both models over the entire GrIS domain and at each time step. Thus, by taking the temporal variations in the simulated time series into account, the spatial RMSE makes it possible to assess the model's performance over both the entire geographical domain and the time period under consideration. It is computed as follows:

$$\begin{array}{}\text{(35)}& \mathrm{\text{SRMSE}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_t\times N}}\sum _{t=\mathrm{1}}^{N_t}\sum _{i=\mathrm{1}}^N{\left(X_{i,\mathrm{OR}}\left(t\right)-X_{i,\mathrm{MAR}}\left(t\right)\right)}^{\mathrm{2}}}.\end{array}$$

X_i,OR(t) and X_i,MAR(t) are respectively the ORCHIDEE-ICE and MAR variables at each grid point i and each time step t. N_t is the number of time steps. In contrast to the RMSE, we used the daily simulated values to compute the SRMSE.

While the RMSE and SRMSE give an indication of the magnitude of the absolute difference between both models, it is also important to calculate the area-weighted average bias (hereafter, areal-mean bias) of each grid point in order to examine whether the model variables simulated by ORCHIDEE-ICE are underestimated (negative bias) or overestimated (positive bias). This bias (MB) is given by

$$\begin{array}{}\text{(36)}& \text{MB}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N{A}_i\left(X_{\mathrm{OR}}\left(i\right)-X_{\mathrm{MAR}}\left(i\right)\right)}{\sum _{i=\mathrm{1}}^N{A}_i}},\end{array}$$

where A_i is the surface area of each grid point.

Finally, we also examined the probability density functions (PDFs) and performed a Cramer–von Mises (CVM) test (Anderson, 1962) to compare the MAR and ORCHIDEE-ICE distributions of a given variable. The CVM test integrates the quadratic differences between the two models over the whole distributions (including the tails of the distributions). In this sense, it is more powerful and more sensitive to departures from the reference distribution (i.e. MAR) than the widely used Kolmogorov–Smirnov test (Stephens, 1970), which is based on the absolute value of the greatest distance between the two distributions.

5 Results

5.1 Evaluation against MAR for standard albedo parameters

Figures 2 to 4 display the spatial distribution of runoff, sublimation, and refreezing simulated by MAR (Figs. 2a, 3a, 4a) and by ORCHIDEE-ICE in the STD-3L (Figs. 2b, 3b, 4b) and STD-12L (Figs. 2c, 3c, 4c) experiments. The main runoff areas simulated with MAR are located on the western edge; however, to some extent, runoff occurs in all peripheral areas of the ice sheet (Fig. 2a). Locations of the ablation zones are well represented in ORCHIDEE-ICE but are limited to a very narrow band, especially in the STD-3L simulation (Fig. 2b). Increasing the number of snow layers favours the inland expansion of the ablation areas on the western and northern margins (Fig. 2c). However, this expansion remains too restricted compared to MAR (Fig. 2a). In the ablation areas, differences in the amount of runoff exceed 1.5 mm d⁻¹ (Fig. S1). Integrated over the whole ice sheet (Table 2), the runoff values computed in the STD-3L (152 Gt yr⁻¹) and STD-12L (205 Gt yr⁻¹) experiments for the 2000–2019 period are respectively 59 % and 45 % lower compared to those of MAR (375 Gt yr⁻¹). As a consequence of the considerably smaller amount of runoff in ORCHIDEE-ICE and thus of surface meltwater, refreezing is also much lower (Table 2) and less extended (Fig. 3a–c) compared to that in MAR. It can be noted, however, that the disagreement is less pronounced with the STD-12L experiment (Fig. S2), which underlines the benefit of increasing the number of snow layers.

Figure 2 Spatial distribution of the runoff (in mm d¹) averaged over the 2000–2019 period and simulated with MAR (a) and the ORCHIDEE-ICE model (b–e) using the 3-layer snow scheme and the standard albedo parameters (b), the 12-layer snow scheme and the standard albedo parameters (c), the 12-layer snow scheme and the albedo parameters optimized using a data assimilation technique (Raoult et al., 2023) and a previous version of the ORCHIDEE-ICE model (d), and the 12-layer snow scheme and the albedo parameters obtained after manual tuning (e).

Table 2 Simulated values of SMB, runoff, sublimation, and refreezing integrated over the entire Greenland ice sheet and averaged over the 2000–2019 period (second column). Evaluation of simulated SMB and SMB components is done with respect to MAR outputs using the values of the root-mean-square error (third column), areal mean bias (fourth column), and spatial root-mean-square error (fifth column). Bold formatting refer to the reference (MAR) and manual tuning (OPT-12L) experiments.

Figure 3 Same as Fig. 2 but for the simulated refreezing (in mm d⁻¹).

Large differences between MAR and ORCHIDEE-ICE also arise regarding sublimation (32 and 33 Gt yr⁻¹ in the STD-3L and STD-12L experiments respectively against 82 Gt yr⁻¹ for the 2000–2019 period in MAR). This feature concerns the entire ice sheet but is even more striking in peripheral areas (Figs. 4 and S3). In central Greenland, differences are smaller, but ORCHIDEE-ICE simulates a little condensation (Fig. 4), whereas MAR does not.

The differences in simulated runoff and in sublimation between MAR and ORCHIDEE-ICE translate into overestimated SMB values simulated with ORCHIDEE-ICE (504 and 450 Gt yr⁻¹ in STD-3L and STD-12L against 286 Gt yr⁻¹ in MAR; see also Figs. 5 and S4). Since inland regions are dominated by the accumulation signal, which is provided by the MAR outputs, the SMB anomalies are primarily driven by differences in the ablation components occurring at the edges of the ice sheet and exceed 2 mm d⁻¹ in most parts of the western and southeastern margins.

An important conclusion that can be drawn from these results is that the use of a better-resolved snow-layering scheme (12-layer as opposed to a 3-layer snow scheme) reduces the mismatch between MAR and ORCHIDEE-ICE. This is mainly illustrated by the integrated SMB and runoff values that are respectively ∼35 % lower and ∼11 % higher in STD-12L, translating into reductions in RMSE values (∼19 % and ∼17 % for SMB and runoff respectively; see Table 2), areal mean bias (∼25 % and ∼24 % respectively), and, to a lesser extent, of the spatial RMSE (∼8 % for both SMB and runoff). Nevertheless, the differences compared with MAR are still too large for the model to be used as a reliable tool to compute the GrIS SMB.

Figure 4 Same as Fig. 2 but for the simulated sublimation (in mm d⁻¹). Negative values indicate condensation.

Figure 5 Spatial distribution of the GrIS SMB simulated with MAR (in mm d⁻¹) and averaged over the 2000–2019 period (a). Differences in the GrIS surface mass balance between MAR and the ORCHIDEE-ICE model (b–e) with the standard parameter values of the albedo parameterization and the three-layer snow scheme (b). Panels (c)–(e) correspond to simulations performed with the updated 12-layer snow scheme for standard values of the albedo parameters (c), optimized values of the albedo parameters (d), and values of the albedo parameters obtained after manual tuning (e).

5.2 SMB and runoff for modified albedo parameters

5.2.1 Impact of optimized albedo parameters

As snow is a highly reflective medium, small changes in albedo may produce large changes in the surface energy balance and thus in the SMB. In the GrIS interior, there is generally quite-good agreement between the summer albedo computed by MAR and the standard ORCHIDEE-ICE simulations (i.e. STD-3L and STD12-L experiments; Figs. 6b and c and S5) with slight negative anomalies of less than 0.05. Negative anomalies ($\sim -\mathrm{0.1}$) also appear, mainly in the northern part of the ice sheet, but with only a few consequences to surface melting owing to the very cold conditions in this region. However, on the western margin where most of the melting takes place, larger snow albedo values are found in ORCHIDEE-ICE. This leads to underestimated surface temperatures compared to MAR (Fig. 7) and, thus, to undervalued runoff that may explain some of the discrepancies between MAR and ORCHIDEE-ICE. There are also differences between the observations provided by MODIS retrievals and the MAR albedo (Fig. 6a and f), especially in the northern and southern parts and on the western margin. On the other hand, the summer albedo computed in the STD-3L and STD-12L experiments (Fig. 6g and h) is generally too low in the interior of the ice sheet and too high on the western margin, with differences from 0.05 to 0.15.

Figure 6 Left – spatial distribution of the summer (JJA) albedo computed with MAR (a) and MODIS (f) and averaged over the 2000–2017 period. Right – differences between the albedo computed with ORCHIDEE-ICE and MAR (b–e) and between ORCHIDEE-ICE and MODIS (g–j) for the 3-layer snow scheme and the standard albedo parameters (b, g), the 12-layer snow scheme and the standard albedo parameters (c, h), the albedo parameters inferred from a data assimilation technique and using a previous version of the ORCHIDEE-ICE model (d, i), and the albedo parameters obtained after manual tuning (e, j).

As mentioned in Sect. 3.2, we investigated the sensitivity of the SMB and its components to the albedo. We first performed an ORCHIDEE-ICE experiment (ASIM-12L) with the optimized albedo parameters inferred from Raoult et al. (2023). Figure 6i illustrates how the representation of the albedo has been improved in the ASIM-12L experiment compared to STD-12L (Figs. 6h, S5 and S8). Model–data discrepancies are now reduced, with differences lower than 0.05 except in the northernmost parts of the ice sheet. The RMSE decreased by ∼26 % (Table 3), which is quite consistent with Raoult et al. (2023). The ablation areas are now better represented (Fig. 2d) due to increased surface temperatures (Fig. 7c) as a result of lower albedo values on the western margin (Fig. 6i).

Table 3 Albedo RMSE values (second column), areal mean biases (third column), and spatial RMSE (fourth column) with respect to MODIS (top) and MAR (bottom).

However, despite the smaller mismatch between modelled ASIM-12L albedo and MODIS retrievals and the better representation of the ablation areas, the simulated amount of runoff (217 Gt yr⁻¹) integrated over the whole GrIS has only slightly been improved with respect to STD-12L (Fig. 2d) and remains quite different from the MAR outputs (Fig. 2a). In addition, the simulated SMB (466 Gt yr⁻¹) has even been slightly degraded (Fig. 5a and d) due to negative temperature anomalies in central Greenland extending to the southern tip (Fig. 7c), resulting from slightly higher albedo values compared to MAR and MODIS (Fig. 6a, i).

Figure 7 Spatial distribution of the snow temperature differences with respect to MAR averaged over the 2000–2019 period (in °C) simulated for the STD-3L (a), STD-12L (b), ASIM-12L (c), and OPT-12L (d) experiments.

The low performance for the SMB computation in ASIM-12L is not solely due to a small amount of runoff but also due to strong negative values of sublimation (i.e. large condensation) over central Greenland (Fig. 3d) resulting in an average level of 5 Gt yr⁻¹ over the entire ice sheet compared to 82 Gt yr⁻¹ in MAR (Table 2). In the ASIM-12L experiment, the albedo in the central GrIS region is slightly higher (up to 0.05) than the albedo retrieved from MODIS (Fig. 6i), while the albedo computed with MAR is slightly lower (Fig. 6a and f). This explains why the ASIM-12L surface temperatures are smaller than those simulated with MAR. This can lead, therefore, to lower saturation pressures that can drop below the dew point and thus produce solid condensation. This result highlights the key influence of the albedo on surface processes and, in particular, illustrates how a small departure from observations may lead to strong biases in sublimation estimates.

5.2.2 Manual tuning

As mentioned in Sect. 3, we have not yet performed a data assimilation experiment to calibrate the new 12-layer ES model, given the computational cost of such an experiment. Instead, we chose to follow a trial-and-error approach. As runoff dominates the SMB signal, our primary objective was to improve the runoff computation by reducing the summer albedo values in the main ablation areas (i.e. the western margin). Given the number of albedo parameters, several options are available to achieve this:

• lowering the albedo of aged snow (A_aged) and/or the albedo of fresh snow (A_aged+B_dec);

• modifying the parameter controlling the decay rate of snow albedo (τ_dec);

• increasing snow age by changing the parameters related to snow ageing, such as the minimum snowfall thickness needed to reset snow age to zero (δ_c), the tuning parameters ω₁, ω₂ (see Eq. 10), and the maximum snow age (τ_max);

• changing the ice albedo (α_ice) because it can also affect SMB and runoff computation if the snowpack melts entirely during summer months in some places and gives rise to bare ice.

Owing to the various influences of the albedo parameters, we had to find a compromise to lower the albedo in ablation areas and improve the computation of runoff and SMB while keeping reasonable albedo values in the GrIS interior. Among the values we tested for each of the parameters, the set of parameters providing the best agreement with MAR outputs (for SMB and SMB components) is highlighted in bold in Table 1 (OPT-12L experiment). Compared to the ASIM-12L experiment (Figs. 6i, S5, S8), the albedo mismatch between ORCHIDEE-ICE (OPT-12L experiment) and MODIS is amplified, especially along the western margin and in the northern sector, with differences reaching 0.25 and 0.3 respectively (Fig. 6j). Nevertheless, these results were expected since our manual tuning was designed to increase the magnitude of the ablation components (especially runoff) and to decrease the SMB and, therefore, to lower the albedo values with a direct impact on surface temperatures, hence surface melting and sublimation.

5.2.3 Impact on SMB components

Using the new set of albedo parameters obtained with the manual tuning approach, the ablation areas are now much more extended than those simulated in the STD-12L experiment (Fig. 2c and e). Compared to MAR (Fig. 2a), they are even wider in the northern part due to increased surface temperatures (Fig. 7d) in response to lower albedo values (up to −0.25). The total amount of runoff averaged over the 2000–2019 period is now 336 Gt yr⁻¹ (against 375 Gt yr⁻¹ in MAR). For the OPT-12L experiment, the RMSE value decreased by ∼40 % compared to STD-12L (Table 2). In the same way, the sublimation (52 Gt yr⁻¹) and refreezing (158 Gt yr⁻¹) match better with MAR (Table 2). In particular, condensation over central Greenland has been considerably reduced, notably with respect to ASIM-12L, but sublimation is still underestimated along the GrIS edges and in the southern part (Fig. 4e). The increase in refreezing (with respect to STD-12L and ASIM-12L) in the GrIS interior (Fig. 3e) is likely linked to lower summer albedo values (Fig. 6e and j) leading to a smaller amount of melting compensated by refreezing. In the main ablation areas, more refreezing is produced and thus there is better agreement with MAR, although it is still insufficient.

These results for the SMB components are evidently associated with an improved representation of the SMB itself (Fig. 5e), which now reaches 301 Gt yr⁻¹ (286 Gt yr⁻¹ obtained with MAR). Indeed, the RMSE and the spatial RMSE values have been reduced by ∼41 % and 10 % respectively for the SMB (∼28 % and 9 % for the runoff) compared to the STD-12L experiment (Table 2). An even more striking result concerns the areal mean bias, which has been lowered by 1 order of magnitude. These improvements are also illustrated in Fig. 8, which displays the monthly mean values for each grid point of the SMB components simulated with ORCHIDEE-ICE as a function of the same MAR variables (see for example the correlation coefficient for both SMB and runoff for the OPT-12L experiment). However, our results are less conclusive for sublimation and refreezing. Although the areal-mean bias and the RMSE values indicate a better match between the OPT-12L and the MAR simulations, the spatial RMSE values are larger compared to the three other ORCHIDEE-ICE experiments, suggesting lower temporal consistency between OPT-12L and MAR. In addition, the correlation coefficients for sublimation and refreezing are also smaller (Fig. 8). On the other hand, the best overlaps between the probability density functions in MAR and the ORCHIDEE-ICE experiments are undoubtedly obtained for OPT-12L, as shown in Figs. S6–S7 and the scores of the CVM tests reported in Table S1.

Figure 8 Representation of the simulated SMB (first row), runoff (second row), sublimation (third row), and refreezing (fourth row) simulated with ORCHIDEE-ICE as a function of the same MAR variables: STD-3L (first column), STD-12L (second column), ASIM-12L (third column), and OPT-12L (fourth column). The different points represent the monthly mean values over the period of 2000–2019 for each of the grid points. The regression line is displayed in red (R is the regression coefficient), and the line y=x is in black.

Despite these encouraging results, it is important to underline that the improved SMB simulation in OPT-12L is achieved through the albedo reduction, and therefore, to some extent, comes from error compensation. However, the reduced albedo also makes it possible to compensate for the effect of some missing mechanisms, such as the lack of consideration of snow–atmosphere interactions or the absence of an explicit representation of snow metamorphism, which has a direct impact on the density profile, the albedo itself, and the temperature profile.

5.3 Vertical temperature and density profiles

To go a step further and gain a better understanding of the above results, it is also important to explore the internal processes of the snowpack. To achieve this, we chose to focus on the vertical temperature and density profiles. Figure 9 shows the snow temperatures simulated by ORCHIDEE-ICE as a function of the MAR snow temperatures at 20cm and 1m depths in the snowpack. These plots show that the temperatures simulated in STD-3L, STD-12L, and OPT-12L behave in approximately the same way as those of MAR. In the first 20 cm, ORCHIDEE-ICE is slightly warmer than MAR for temperatures between −30 and −10 °C, despite a few slightly colder grid points appearing in the range of −20 to −10 °C. The ASIM-12L experiment presents the best agreement with MAR, although it has slightly lower temperatures. These features directly reflect the behaviour of surface temperatures (Fig. 7) that strongly influence the upper snowpack layers. Another key point arising from these plots is the very good agreement between MAR and ORCHIDEE-ICE for temperatures above −10 °C. This suggests that the potential runoff that could occur in the first tens of centimetres of the snowpack should not be affected so much.

However, the departure from MAR increases with snow depth, especially for low temperatures. For example, at 1 m depth, differences of 3–4 °C are obtained (Fig. 9) and may exceed 5 °C for deeper levels (not shown). These enhanced differences from MAR are likely due to a positive feedback related to the thermal conductivity (see Eq. 33): as snow temperature increases by 1 °C in a given layer, the thermal conductivity increases by 1 order of magnitude.

As pointed out by Domine et al. (2019), the snow thermal regime and snow density are strongly coupled. As an example, they mentioned the work of Fréville (2015), who showed that an error of 1 °C in the surface temperature can lead to errors in snow density of 100 kg m⁻³. Our experiments show that for a depth of 20 cm, the higher the surface temperature, the lower the snow density on average (Fig. 10). On the other hand, in the ASIM-12L experiment, snow temperatures are lower compared to the three other ORCHIDEE-ICE experiments, and snow densities are larger. This contradicts a number of studies (e.g. Kojima, 1967; Anderson, 1976; Mizukami and Perica, 2008), which have shown that in a warmer snowpack, snow grains become rounded and are more prone to be compacted more easily, hence leading to an increase in snow density. However, in our model this process cannot be reproduced, as snow metamorphism is only accounted for through snow ageing. Conversely, in deeper layers, the model is more effective at densifying (Fig. 10), in line with the fact that warmer snow becomes more plastic and compacts more easily. In particular, between 20 cm and 1 m in depth, the RMSE computed between OPT-12L and MAR has been reduced from 79.63 to 30.22 kg m⁻³. Beyond 600 kg m⁻³, the ORCHIDEE-ICE densities are generally below those of MAR because the maximum density is fixed to 750 kg m⁻³ (see Sect. 2). However, the comparison of our results on snow density with those of MAR should be viewed with caution because, to the best of our knowledge, the snow density simulated by MAR has not been evaluated against available observations.

Figure 9 Representation of the ORCHIDEE-ICE simulated snow temperatures at 50 cm (left) and 1 m in depth (right) as a function of the MAR snow temperatures. The different points represent the monthly mean values over the period of 2000–2019 for each grid point. The regression line is displayed in red (R is the regression coefficient), and the line y=x is in black.

Figure 10 Same as Fig. 9 but for snow density.

5.4 SMB evolution: impact of the climate forcing

The results presented in the previous sections were averaged over the 2000–2019 period (for SMB and the SMB components) and over the 2000–2017 period (for the albedo). In this section, we present the temporal evolution of the SMB between the years 2000 and 2019 (Fig. 11). Figure 11 shows that regardless of the ORCHIDEE-ICE experiment under consideration, the evolution of the yearly integrated SMB is in accordance with the evolution simulated by the MAR model. In particular, the years in which extreme melting events were recorded (such as 2012 and 2019) are perfectly well represented (Bennartz et al., 2013; Tedesco and Fettweis, 2020). As expected, the best agreement with MAR is obtained for the OPT-12L experiment as a result of the calibration of the albedo parameters.

When forced by the ERA-5 meteorological fields and using the manually tuned parameters, ORCHIDEE-ICE simulates higher SMB values and a lower runoff (Fig. 11 and Table 2), especially during the first period of the time series (2000–2008). However, the evolution of the yearly integrated SMB in the ERA5-12L experiment follows exactly the same interannual variations as for the OPT-12L experiment forced with MAR (Fig. 11). This indicates that the surface climate simulated by MAR is close to that derived from the ERA-5 products. Moreover, in a comparative study of the ERA-5 reanalyses, Arctic system reanalysis, and MAR performance, Delhasse et al. (2020) showed that MAR outperforms ERA-5 for near-surface temperatures when compared to observations from automatic weather stations. As the surface melt and thus the SMB largely depend on near-surface temperatures, there is, therefore, a strong interest in using MAR to force our snow model and to compare its performance to that of MAR.

Figure 11 Evolution of the yearly surface mass balance of the Greenland ice sheet simulated with MAR (orange), ORCHIDEE-ICE forced by MAR outputs (STD-3L and STD-12L – blue, solid, and dashed lines; OPT-12L – purple line), and ORCHIDEE-ICE forced by the ERA-5 reanalyses (green line).

In this paper, we have so far limited the comparison of our results to those of MAR. However, as mentioned in Sect. 4, we also evaluated the simulated SMB with 353 daily SMB observations available from the PROMICE database over the 2000–2019 period (Machguth et al., 2016; Mankoff et al., 2021a, b). In addition, it is also interesting to evaluate our model results against observations when ORCHIDEE-ICE is forced by climatic fields independent from MAR outputs. To address this issue, we plotted the modelled SMB for OPT-12L, ERA5-12L, and MAR for the grid points located closest to the observation sites as a function of the PROMICE SMB measurements (Fig. 12). We also provided statistical elements for the comparison between MAR, the five ORCHIDEE-ICE experiments, and the SMB observations (Table 4). This model–data comparison confirms the conclusions we reached when evaluating the performance of our model against MAR outputs, namely the significant improvement in our results when moving from STD-3L to OPT-12L. Moreover, although the bias between the OPT-12L SMB and the observed SMB is twice as high as for MAR, the model–data correlation is of the same order of magnitude as for MAR (Table 4).

Figure 12 Simulated SMB in the OPT-12L experiment and in MAR as a function of the observed SMB from the PROMICE network. As the observed SMB values are not all available over the same time interval, the measurements are given in metres of water equivalent (mWE); 353 observations were available over the 2000–2019 period. Each simulated SMB value corresponds to the grid points located closest to the observation sites. The red line is the regression line, with R being the correlation coefficient, and the black line indicates the line y=x.

The ERA5-12L experiment also produces good agreement with the observations. Despite a lower correlation coefficient than for MAR and OPT-12L, the mean bias is on the same order of magnitude as that of MAR, and the RMSE for the SMB obtained is the lowest for any of the experiments. It is clear that the SMB simulated in the experiments forced by MAR is partly driven by the climate simulated by MAR itself (for the accumulation component). However, the results obtained with ERA5-12L clearly show that the behaviour of our model is consistent regardless of the climate forcing used. Nevertheless, it should be noted that the resolution of ORCHIDEE-ICE corresponds to that of the model used as a forcing. For ERA5-12L, the resolution is about twice as fine as for the experiments forced by MAR (20 km × 20 km). Thus, to make our comparison between ERA5-12L, MAR, and/or OPT-12L more robust, we should have used MAR with a resolution of 10 km × 10 km. We cannot therefore rule out that the results for OPT-12L would then have provided a better comparison with the PROMICE data than ERA5-12L.

Table 4 Comparison of the simulated SMB in MAR, STD-3L, STD-12L, ASIM-12L, and OPT-12L with the SMB observations from the PROMICE network. The bias is computed as the average between modelled and observed SMB for each grid point. Note that the values of the bias and the RMSE are given in mWE, as the observed SMB values are not all available over the same time interval.

6 Discussion and concluding remarks

The land surface component of the IPSL ESM used for CMIP6 included a three-layer snowpack model operating over continental surfaces. However, this snow scheme was not adapted to glaciated surfaces, which is a major drawback and makes it impossible to compute the surface mass balance over ice sheets or glaciers. The aim of this paper was therefore to present the new developments made to adapt the snow model to ice-covered areas and to document its performance. Our first step was to calibrate the snow albedo parameterization over the Greenland ice sheet. In order to have a set of climate variables covering the whole ice sheet, we chose to force our model by the atmospheric outputs of the MAR regional model, which shows very good performance at simulating the surface climate and thus offers undeniable advantages for the representation of the physical processes related to snow and ice, in particular surface melting (Delhasse et al., 2020). We have shown that the ablation-related processes are highly dependent on the choice of the albedo parameters. The set of parameters obtained after manual tuning (OPT-12L experiment) provides good agreement between the SMB computed in ORCHIDEE-ICE and MAR. However, as outlined in Sect. 5.2.3, this improvement is mainly the result of albedo lowering. The summer albedo computed with this set of parameters has been degraded compared to MAR and MODIS and to the albedo computed in the ASIM-12L experiment (based on the MODIS-optimized albedo parameters), as shown in Table 3 and in Figs. 6i–j, S5, and S8. While the RMSEs computed between ORCHIDEE-ICE and MAR for SMB and runoff have been reduced by ∼39 % and ∼33 % respectively from ASIM-12L to OPT-12L, the RMSE for albedo has increased by 47 % (Table 3). The mismatch between MODIS retrievals and the OPT-12L albedo is mainly observed in the northernmost part of the ice sheet and, to a lesser extent, on the western edge.

A more objective method would have been to perform a data assimilation experiment similar to the one presented in Raoult et al. (2023) using the new version of the ORCHIDEE-ICE model. However, albedo is not the only important parameter governing the snowpack evolution. The albedo parameters inferred from the Raoult et al. (2023) optimization greatly improve the representation of the albedo but degrade the other model outputs compared to those obtained with the manually tuned albedo parameters. This is most likely because their optimization overfits the albedo retrievals without applying constraints to other processes, strongly impacting the SMB components and controlling the state of the snowpack (e.g. snow compaction, snow density, snow viscosity). This supports the recommendation for a multi-objective optimization using not only albedo data but also vertical temperature and density profiles, as well as SMB observations. Since this type of approach is highly time-consuming, it has not yet been undertaken but could be the objective of a future study.

However, the reduction in albedo in the current ORCHIDEE-ICE version can compensate for missing processes. For example, snow drift, transmission of solar radiation, or the effect of light-absorbing particles on the albedo are ignored. Metamorphism is not explicitly represented, although its effect on the albedo and the vertical density profile are accounted for (albeit in a crude manner) through the snow ageing function f_age (Eq. 7) and the ψ_snow function (Eq. 17) respectively.

In the GrIS, snow erosion has often been considered a second-order component of mass loss in ablation areas compared to meltwater. However, in the ice-sheet interior, sublimation and snow erosion are dominant processes in removing mass from the surface and may have, therefore, a significant impact on SMB (van Angelen et al., 2011).

Taking into account the transmission of solar radiation within the snowpack can lead to a warming of the internal layers, with higher temperatures near the surface and lower temperatures at depth due to the exponential decrease in heat transfer. This results in a temperature gradient that influences the metamorphism of snow grains and thus accelerates densification (Colbeck, 1983). We showed that the ORCHIDEE-ICE temperatures inside the snowpack were higher than those simulated by the MAR model. A likely hypothesis to explain this behaviour relies on the reduction in albedo, which leads to excessively high surface temperatures. However, it is important to note that heat transfer can promote snow melting, which in turn can percolate at depth and refreeze, affecting both the runoff and the vertical structure of the snowpack through changes in density (Colbeck, 1983). Quantifying all these processes requires, therefore, the proper representation of solar absorption, which is itself strongly dependent on snow optical properties (Warren, 1982) and, therefore, on snow grain size (Libois et al., 2013). Since metamorphism is not explicitly represented in the model, we assumed that representing solar absorption was not a priority in our modelling approach, even if this choice is debatable. However, in the near future, a more sophisticated albedo scheme based on a radiative-transfer model accounting for light-absorbing particles and snow grain size (Kokhanovsky and Zege, 2004) will be implemented in the ORCHIDEE-ICE model. This will allow us to represent the backward- and forward-scattering processes as well as light absorption.

In addition, there are also structural deficiencies related to the fact that in ORCHIDEE-ICE, a single energy balance is computed in one grid cell. This is detrimental for the albedo computation, especially at the edges of the ice sheet where several surface types may coexist in a 20 km × 20 km mesh. However, the implementation of a multi-tile energy balance is currently under development.

Finally, as our simulations have been run in offline mode, the snow feedback onto the atmosphere has not been taken into account, in contrast to the MAR model that is fully coupled to a snow scheme derived from CROCUS (Brun et al., 1989, 1992). Ignoring snow–atmosphere feedback may potentially lead to biases related to surface processes and to an improper representation of the energy and humidity flux exchanges at the snow–atmosphere interface. For example, forcing our model with the atmospheric temperature at 2 m derived from the fully coupled MAR simulation could lead to an underestimation of the energy available at the snow–atmosphere interface, resulting in less snowmelt compared to what is simulated in coupled mode. However, our manual tuning approach aims to limit the potential underestimation of the surface meltwater production. Conversely, any potential bias in the MAR forcing may also affect our results (Dietrich et al., 2024). To overcome this problem, it would have been interesting to force ORCHIDEE-ICE with meteorological fields recorded at the automatic weather stations. This has not been done in this study because the meteorological fields required to force ORCHIDEE-ICE were not all available at the PROMICE stations and because our first objective was to obtain a reasonable estimate of the SMB and its components at the scale of the entire GrIS.

Despite the potential improvements that could still be made to ORCHIDEE-ICE to enhance the model's performance, the developments presented in this paper represent a major step forward. Indeed, they now allow the ice-sheet surfaces to be handled by the land surface model, consistently with all the other surface types, and not by the atmospheric component of the IPSL model (LMDZ) as was the case up to now. In addition, the new snow model can now be applied to the continental glaciers, replacing the very crude snow scheme used previously. Our developments enable us to provide a reasonable estimate of the surface mass balance of the Greenland ice sheet, in very good agreement with that simulated by the MAR model, which was used as a reference in this study. These developments constitute a first step towards the full coupling between the IPSL global climate model and ice-sheet models.

Appendix A

Table A1 List of variables used in ORCHIDEE-ICE and related to snowpack and ice processes.