# Dry snow initialization and densification over the Greenland and Antarctic ice sheets in the ORCHIDEE land surface model

Philippe Conesa, Cécile Agosta, Sylvie Charbit, Christophe Dumas, Simon Beylat, Nina Raoult

► To cite this version:

Philippe Conesa, Cécile Agosta, Sylvie Charbit, Christophe Dumas, Simon Beylat, et al. Dry snow initialization and densification over the Greenland and Antarctic ice sheets in the ORCHIDEE land surface model. 2026. hal-05529885

HAL Id: hal-05529885

https://hal.science/hal-05529885v1

Preprint submitted on 27 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Dry snow initialization and densification over the Greenland and Antarctic ice sheets in the ORCHIDEE land surface model

Philippe Conesa<fnref n="1" />,<fnref n="2" />, Cécile Agosta<fnref n="1" />,<fnref n="2" />, Sylvie Charbit<fnref n="1" />,<fnref n="2" />, Christophe Dumas<fnref n="1" />,<fnref n="2" />, Simon Beylat<fnref n="1" />,<fnref n="2" />, and Nina Raoult<fnref n="3" />

**Correspondence:** Philippe Conesa (philippe.conesa@lsce.ipsl.fr)

**Abstract.** Accurate modeling of the snowpack over ice sheets is essential for quantifying their surface mass balance contribution and their resulting impact on sea level rise. The snowpack evolution is largely governed by surface climate but also by internal processes such as densification. In land surface models, such processes are often represented using formulations developed for seasonal snow, limiting their realism in polar environments. To improve the representation of polar snow in the land surface model ORCHIDEE (Organizing Carbon and Hydrology in Dynamic Ecosystems), the surface component of the IPSL-CM climate model, we implement a series of developments over Greenland and Antarctica. These developments include (1) a new snowpack initialization procedure that generates deep, physically consistent density profiles based solely on latitude and elevation aimed to be included in other snowpack models, (2) a wind-based surface snow density parameterization applicable to both Greenland and Antarctica, and (3) a recalibrated dry-snow densification scheme for snowpack compaction, using observations and 1D and 2D offline simulations. These developments improve the simulation of surface snow density as well as the internal snowpack structure, including both density and temperature, with good agreement with dry-snow observations. Discrepancies still persist in representing the Greenland ablation zones, suggesting that further improvements in surface energy balance processes are needed, with specific attention to snow albedo.

## 1 Introduction

The evaluation of the ice sheets mass balance is crucial to understand their evolution and their potential contribution to sea-level rise (Alley et al., 2005; Box et al., 2022; Otosaka et al., 2023). In Antarctica, ice sheet mass changes are primarily driven by basal melting and ice discharge into the ocean (Rignot et al., 2019), whereas in Greenland, mass loss is dominated by surface processes (van den Broeke et al., 2016; Mouginot et al., 2019). Due to higher atmospheric temperatures and the subsequent increase in surface melt, the mass loss of the Greenland ice sheet (GrIS) has accelerated over recent decades, contributing to approximately 20% of the observed global sea-level rise (Forsberg et al., 2017; Otosaka et al., 2023). In the Antarctic ice sheet (AIS), surface ablation is limited to specific areas mainly located in the Antarctic Peninsula and the ice shelves (Shepherd et al., 2018; Leduc-Leballeur et al., 2020). Overall, the AIS contributes to 10% of the observed global mean sea-level rise since

<footnote n="1">¹Laboratoire des Sciences du Climat et de l'Environnement CEA-CNRS-UVSQ UMR 8212, Paris-Saclay University, CEA Saclay, 91191 Gif-sur-Yvette cedex, France</footnote>

<footnote n="2">²Institut Pierre Simon Laplace, Sorbonne University, 4 place Jussieu, 75252 Paris cedex, France</footnote>

<footnote n="3">³European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading, RG2 9AX, United Kingdom</footnote>
<figure><img src="image_3.png" /></figure>

1993 (Oppenheimer et al., 2022). Under continued surface warming, the mass balance of the ice sheets is projected to further decrease, driven by enhanced surface melt, ice discharges and basal melting (Oppenheimer et al., 2022).

Surface ablation is strongly related to the processes occurring on the surface and within the snowpack overlying the ice sheet. The snowpack is a porous material composed of snow crystals, air and, potentially, liquid water. This porosity enables the retention and refreezing of surface meltwater or rainwater that percolates within the snow layers (Vandcrux et al., 2020; Amory et al., 2024). Over the GrIS, the increased surface melt and the expansion of the melt area (Mouginot et al., 2019; Mousavi et al., 2021) affect the snowpack internal properties, notably inducing an increase in ice concentration, (de la Peña et al., 2015; Vandcrux et al., 2018), a rise in the snowpack temperature driven by an increase in meltwater percolation and refreezing (Polashenski et al., 2014; van den Broeke et al., 2016) and a decrease in pore spaces reducing the retention and refreezing of water (van Angelen et al., 2013; Vandcrux et al., 2019). These changes are closely linked to variations in snow density, a fundamental property that constrains the pore space available for meltwater retention and refreezing. As a consequence, a detailed understanding and modeling of physical processes governing densification is crucial to constrain the current and future mass changes of both ice sheets.

Different stages of dry snow densification can be distinguished, starting with the surface densification that occurs at the top of the snowpack. According to Arctic field studies, the surface density in polar regions is related to the combined effects of wind compaction and the upward transfer of water vapor induced by the high temperature gradient within polar snowpacks (Barrere et al., 2017; Domine et al., 2016). Moreover, frequent strong winds cause large amounts of snow to be transported, compacted, and sublimated, contributing to the denser surface layers (Liston and Sturm, 2002; Sturm et al., 2001). Recent wind tunnel experiments further quantified surface snow densification, showing that it strongly depends on wind speed and the duration of wind exposure (Walter et al., 2024). Similarly, observations during a drifting-snow event in Adélie Land revealed rapid post-depositional snow densification, with density increasing from approximately 200 kg m⁻³ to around 350 kg m⁻³ in a single day (Amory et al., 2021). However, there is a lack of surface observations in polar regions limiting our ability to determine which processes are involved in surface densification.

Below the surface, two main dry-snow densification stages can be distinguished. The first one occurs until density reaches around 550 kg m⁻³ and is characterized by rapid compaction associated with grain settling and packing (Herron and Jr, 1980; Arnaud et al., 1998). Between 550 and 830 kg m⁻³, the densification rate declines and deformation is primarily due to viscous creep and sintering under an increasing overburden pressure (Arnaud et al., 1998). Around 830 kg m⁻³, known as the bubble close-off density, ice forms and air bubbles become trapped. Densification slows significantly as air can no longer escape, and compaction occurs by compressing the trapped air bubbles (Herron and Jr, 1980; Goujon et al., 2003; Stevens et al., 2020). Temperature and accumulation regulate the densification rate by controlling snow viscosity and mass loading, while metamorphism modifies the snow microstructure rather than driving densification itself (Arnaud et al., 1998; Goujon et al., 2003). As a result, densification can take years to centuries under cold or low accumulation conditions. This way, the interior of the ice sheets densifies slowly, whereas warmer and/or higher-accumulation regions (such as the southeastern part of the GrIS or the West AIS) densify more rapidly (van den Broeke, 2008; Ligtenberg et al., 2011; Vandcrux et al., 2020; van Kampenhout et al., 2017).
<figure><img src="image_4.png" /></figure>

Several parameterizations have been proposed to describe snow densification. Early approaches expressed the densification rate as an Arrhenius-type function of temperature, with an additional dependence on accumulation (Herron and Jr, 1980; Barnola et al., 1991). Such parameterizations have been later refined by introducing distinct activation energies for grain growth and vapor diffusion (Arthern et al., 2010; Ligtenberg et al., 2011). This framework has been adopted in polar atmospheric models, including RACMO (Ligtenberg et al., 2011; Kuipers Munneke et al., 2008). Alternative formulations link densification directly to snow microstructural properties, such as grain size and shape (Alley, 1987; Arnaud et al., 1998; Goujon et al., 2003). Other models represent the change in snow density as a function of snow viscosity, inspired by the formulation proposed by Anderson (1976), later implemented and adapted in ISBA-ES (Decharme et al., 2016), ORCHIDEE (Wang et al., 2013; Charbit et al., 2024), CROCUS (Vionnet et al., 2012), and CESM (van Kampenhout et al., 2017). In these models, viscosity is expressed as a function of density and temperature. In more advanced schemes, viscosity also depends on the liquid water content (Decharme et al., 2016) and microstructure properties such as grain diameter and shape (Bartelt and Lehning, 2002; Steger et al., 2017; Vionnet et al., 2012).

60

65

Since many densification rate schemes rely on the physical constraint imposed by the snow mass, their performance is sensitive to how the snowpack thickness is initialized. In Greenland and Antarctica, the ice sheet snowpack can reach several decades of meters (Amory et al., 2024). To form a thick snowpack over polar regions, many models such as ORCHIDEE, CESM, SNOWPACK and IMAU-FDM therefore rely on a spin-up procedure where the snowpack is built through natural accumulation from a snow-free state. Under the assumption of climatic stationarity, this process can take decades to over a century, depending on the local accumulation rate (Charbit et al., 2024; van Kampenhout et al., 2017; Ligtenberg et al., 2011; Keenan et al., 2021; Veldhuijsen et al., 2023). While widely used, such approaches are time-consuming and highlight the need for standardized, physically consistent, and efficient methods to better constrain the initial state of the ice sheet snowpack.

70

75

In this study, we evaluate and improve the capacity of the ORCHIDEE (Organising Carbon and Hydrology in Dynamic Ecosystems) land surface model (Krinner et al., 2005) of the IPSL (Institut Pierre Simon Laplace) Earth System Model (Boucher et al., 2020; Cheruy et al., 2020), to represent snow densification processes over GrIS and AIS. To this end, in Section 3, we develop and test (1) a new automatic initialization method for Antarctica and Greenland snowpacks, generalizable to every snowpack model, (2) a new surface snow density parameterization, common for both ice sheets, and finally (3) a recalibration of the ORCHIDEE densification scheme using the History Matching (HM) method (Williamson et al., 2013), applied to 1D simulations at two reference sites. In Section 4, we then evaluate the calibrated and updated model in 2D offline configurations using observed density profiles across both ice sheets. Finally, we study the impact of these developments on the snowpack and its internal characteristics, such as snow temperature, liquid water retention and refreezing (Section 5).

80

85

# 2 Data and modeling setup

## 2.1 Observations

To improve the representation of snow densification, we use observations of snow temperature, accumulation and density in Greenland and Antarctica obtained from multiple sources, including the SUMup collaborative database (Vandcrux et al.,

90
<figure><img src="image_5.png" /></figure>

2023), the GLACIOCLIM-SAMBA dataset (Favier et al., 2013; Wang et al., 2016) and additional density sources listed in Table 1. For Greenland, we rely only on temperature and accumulation observations recorded after the year 2000 to ensure that our analysis reflects recent climatic conditions, given the significant impact of contemporary climate warming on the GrIS. Temporal restrictions were not imposed on the AIS dataset, as we assume that, at the ice sheet scale, the surface mass balance (SMB) has experienced only minor variations over recent decades (Shepherd et al., 2018).

To compare point observations and model gridded outputs, the observations are averaged on the ORCHIDEE stereographic grid (15 km resolution in Greenland, 35 km in Antarctica). From the SUMup database, we use 3074 observations of snow temperature collected at different locations in Greenland, between 2000 and 2023, and 905 in various sites over Antarctica between 1957 to 2021 (Figure 1a, d). After averaging temperature observations onto the model grid, we obtain 61 points in Greenland and 782 in Antarctica.

The observations of snow accumulation are sourced from the GLACIOCLIM-SAMBA dataset, detailed in Favier et al. (2013) and updated by Wang et al. (2016), with 3625 observations over Antarctica, recorded between 1930 and 2009, and from the SUMup dataset for Greenland with 224 accumulation observations, collected between 2000 and 2022. Over Greenland, we only retain data for which the duration of observations is specified (ranging from a few days to several years). We then compute the grid-averaged accumulation, weighted by the corresponding duration of observations. In total, 61 grid-cell averages are obtained this way. For Antarctica, no grid averaging is applied. To investigate surface snow density, we use density profiles containing measurements within the top 20 cm (560 observations in Greenland and 139 in Antarctica). By averaging the observations on the model grid, we obtain 149 grid cell averages for Greenland and 79 for Antarctica (Fig. 1b, e).

The representation of the vertical densification is assessed using 126 profiles extracted from the SUMup database (79 in Greenland and 49 in Antarctica) and two additional profiles located at Dome C in Antarctica (Leduc-Leballeur et al., 2015) (Fig. 1b, e). Among these observations, 102 reach depths of at least 10 m (79 in Greenland and 23 in Antarctica). We select two study sites characterized by dry snow conditions: Summit (72.58 °N, 38.50 °W; 3254 m a.s.l.) and Dome C (75.06 °S, 123.21 °E; 3233 m a.s.l.). These locations present different climatic conditions. Dome C is characterized by an extremely cold and arid climate with mean annual surface temperatures of -55°C (Buizert et al., 2021) and accumulation rates close to 30 kg m⁻² yr⁻¹ (28.6 kg m⁻² yr⁻¹ observed between 2006 and 2013; Genthon et al., 2016). Summit exhibits a less cold environment, with a mean annual surface temperature observed of -31.4 °C observed between 1998 and 2017 (Vandcrux et al., 2020) and significantly higher accumulation rates of approximately 242 kg m⁻² yr⁻¹ recorded at Summit Camp between 2017 and 2020 (Howat, 2022). At Summit, five density profiles are selected, reaching at least 10 m depth, collected between 2006 and 2017 (Vandcrux et al., 2023). At Dome C, we use two density profiles collected in the 2012-2013 austral summer (Leduc-Leballeur et al., 2015). Figures 1c, 1f represent the observed density profiles at Summit and Dome C respectively. To filter out density measurements fluctuations, we have applied a third-degree polynomial interpolation to the raw density profiles (Fig. 1c,f).
<figure><img src="image_6.png" /><figcaption>**Table 1.** Additional density observations supplementing the SUMup database</figcaption></figure>

<table><thead><tr><td>Author</td><td>Name</td><td>Lat</td><td>Lon</td><td>Elevation [m]</td></tr></thead><tbody><tr><td>Kipfstuhl et al. (2009)</td><td>Kohnen</td><td>-75</td><td>90.07</td><td>2892</td></tr><tr><td>Marquetto et al. (2020)</td><td>West Antarctica</td><td>-79.55</td><td>94.21</td><td>2122</td></tr><tr><td>Ekaykin et al. (2021)</td><td>Ridge B</td><td>-79.02</td><td>93.69</td><td>3400</td></tr><tr><td>Ekaykin et al. (2023)</td><td>Vostok</td><td>-78.27</td><td>106.5</td><td>3480</td></tr><tr><td>Leduc-Leballeur et al. (2015)</td><td>Dome C</td><td>-75.06</td><td>123.2</td><td>3233</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 1.** Mean annual surface temperature for Greenland (a) and Antarctica (d), observed at 10 m depth in the snowpack (circles). Color-filled areas are obtained with the initialization procedure (see Section 3.1). Mean surface density for Greenland (b) and Antarctica (e) for the observed top 20 cm (circles) and computed from the initialization procedure (color-filled areas). Star-shaped markers denote density profiles reaching 10 m depth. Red dashed circles indicate the location of the Summit station ('S') in Greenland (b) and Dome C ('D') in</figcaption></figure>
<figure><img src="image_8.png" /></figure>

Antarctica (e). Observed density profiles at the Summit (c) and Dome C (f) station. Grey lines correspond to density observations and colored curves represent the fitted profiles obtained from the observed density profiles. The dashed pink curve shows the density profile computed by ORCHIDEE (OR OLD) with the standard version of the model (Section 2.3).

## 2.2 The Explicit Snow module in the ORCHIDEE surface model over ice sheets

125

The ORCHIDEE model simulates the water, carbon and energy exchanges between the land surfaces and the atmosphere. ORCHIDEE includes the Explicit Snow module developed to represent the snowpack over continental surfaces (Wang et al., 2013) and recently adapted for glaciers and ice sheets in order to simulate their surface mass balance (Charbit et al., 2024). This module represents the snowpack as a one-dimensional system with 12 layers of snow overlying 8 layers of ice. It simulates surface and internal snow physical processes, including snow densification, heat diffusion, snow aging, surface melting, sublimation, liquid water percolation and refreezing. The snow mass evolves as a function of the liquid and solid precipitation rate, sublimation and the runoff of liquid water that cannot be retained in the snowpack. The surface energy balance is computed using the albedo scheme proposed by Chalita and Le Treut (1994). The ice layers beneath the snowpack are considered as an infinite reservoir, which is only used to compute the contribution of ice melting to runoff in case of complete snow removal during the ablation season. The snow and ice modules have been extensively described in Charbit et al. (2024). Both modules were implemented in ORCHIDEE version 2 (Cheruy et al., 2020). Here, we rely on the Explicit Snow scheme from Charbit et al. (2024) updated in ORCHIDEE version 4. Changes between both ORCHIDEE versions do not affect the representation of snow processes over ice sheets. Prior to the new developments presented in this article, we have addressed several limitations identified in the previous version of Explicit Snow. Previously, refreezing was computed before percolation of liquid water and was therefore limited to the meltwater produced and retained within a given layer. We have revised this approach to allow refreezing to occur after percolation and to depend on the amount of energy available to refreeze, following the methodology implemented in CROCUS (Vionnet et al., 2012). These modifications enable the snowpack to refreeze as much liquid water as its energy balance allows, rather than being restricted by the maximum water retention capacity of individual layers. As a result, the snowpack can now refreeze a larger amount of liquid water as long as sufficient energy is available. Additionally, we increased the maximum snow density in the model from 750 kg m⁻³ to 917 kg m⁻³ to represent the transition from snow to ice densities. Finally, we modified the maximum snow mass threshold from 3000 kg m⁻² to 10000 kg m⁻². This adjustment ensures the formation of a deeper snowpack to better represent percolation, refreezing and runoff.

130

135

140

145

## 2.3 Experimental setup and forcing

We use the ORCHIDEE model in offline mode with a 30 minute time step. The atmospheric forcing is provided by the polar-oriented regional atmospheric model MAR (Fettweis et al., 2017; Kittel et al., 2021), coupled to the SISVAT (Soil Ice Snow Vegetation Atmosphere Transfer) land surface scheme (Ridder and Schayes, 1997), that includes a multi-layered snow model based on the CROCUS model (Brun et al., 1989; Agosta et al., 2019). We set up two configurations, one for GrIS and one for the AIS using the same stereographic grid as in the atmospheric forcing. For Greenland, we use MAR version 3.13 (Fettweis et al., 2023) with a horizontal resolution of 15 km, covering 1980–2019. For Antarctica, MAR version 3.12 is employed, spanning 1980-2019 with a 35 km horizontal resolution (Lambin et al., 2022). Since the MARv3.13 forcings were

150

155
<figure><img src="image_9.png" /></figure>

not available for Antarctica at the time this study was initiated, we use MAR version 3.12 for Antarctica, spanning 1980-2019
with a 35 km horizontal resolution (Lambin et al., 2022). These MAR simulations were forced every 6 hours at their lateral
boundaries by the ERA5 reanalysis (Hersbach et al., 2020) for the period 1980-2019. Considering the albedo representation in
ORCHIDEE, we use for the GrIS the assimilated albedo parameters from MODIS satellite albedo observations (Raoult et al.,
2023), referred also as the ASIM12L parameters in Charbit et al. (2024). Additional simulations, presented in Section 5.2, were
performed with manually tuned albedo parameters, referred to as OPT12L in Charbit et al. (2024). Because albedo is more
spatially uniform and temporally stable over the AIS than over the GrIS, we manually adjusted the albedo parameters for the
AIS, to accurately represent the surface conditions at Dome C. All AIS simulations use the same albedo values (Table S1).
The reference simulation, presented as OR OLD, uses the snow internal parameters introduced in Charbit et al. (2024). This
configuration shows an underestimation of density at the surface and an excessive compaction at depth (Fig. 1).

## 2.4 History Matching calibration method

History Matching (HM) is a method for calibrating parameters (Couvreux et al., 2021; Hourdin et al., 2021), implemented within the ORCHIDEE Data Assimilation Systems (ORCHIDAS MacBean et al., 2022; Raoult et al., 2024). This procedure is based on Gaussian process emulation and is designed to constrain the free parameters of complex models. It has been applied to calibrate atmospheric, ocean and Earth system models (Couvreux et al., 2021; Hourdin et al., 2021; Villefranque et al., 2021; Williamson et al., 2017; Hourdin et al., 2023), as well as ice sheet models (McNeall et al., 2013) and land surface models (Baker et al., 2022; McNeall et al., 2024). Originally developed to calibrate model parameters at the process level, HM operates in ORCHIDAS using offline ORCHIDEE simulations in a single-column model framework. This approach has previously been applied to the ORCHIDEE model to simulate net ecosystem exchange and latent heat fluxes (Raoult et al., 2024). This calibration approach differs from the traditional data assimilation method, which primarily focuses on identifying the optimal parameters that best match observations by minimizing a cost function. In contrast, HM aims to identify regions of the parameter space that are "not plausible" given the constraints predefined through observations, i.e. leading to model outputs inconsistent with observations. The calibration starts by sampling the prior parameter space using a Latin Hypercube design (Morris and Mitchell, 1995), to efficiently explore the space and maximize the information gained from each parameter combination. Following the recommendation of Loeppky et al. (2009), we use a sample size equal to $N = 10 \cdot P$, where $P$ is the number of parameters being calibrated. One-dimensional ORCHIDEE simulations are performed for each sampled set of parameters. Iteratively, Gaussian process emulators are then constructed to capture the relationship between model parameters and scalar metrics, specifically designed to characterize the physical processes under study. The construction of the Gaussian process emulators is done following the methodology presented in Williamson et al. (2013) and the R-scripts from Couvreux et al. (2021). Each metric is associated with its own emulator, which is then used to predict metric values across the full parameter space. This enables the identification of parameter regions that are inconsistent with observations using an implausibility criterion *I*, used to reduce the parameter space. The implausibility is defined as follows:
<figure><img src="image_10.png" /></figure>

$$I(x) = \frac{|m - E[G(x)]|}{\sqrt{\sigma_{obs}^2 + \sigma_{OR}^2 + Var[G(x)]}} \quad (1)$$

with $x$ being a sample of the parameter space, $m$ the reference observed metric, $E[G(x)]$ the value predicted by the emulator $G(x)$. The implausibility accounts for the aggregated uncertainty, computed as the square root of the sum of squared contributions from three sources: the observational ($\sigma_{obs}$), model ($\sigma_{OR}$), and emulator ($Var[G(x)]$) uncertainties. This criterion is computed independently for each metric, and used to define the subregion of the parameter space that is not inconsistent with the reference observational metrics, also called the NROY space for "Not Ruled Out Yet". The NROY space corresponds to the region of the parameter space where the implausibility is smaller than a defined threshold. For the first iterations (waves) of the algorithm, it is common to set the implausibility cut-off to 3, following the 3$\sigma$ rule stating that 99.7% of the probability density for any unimodal distribution is comprised within 3 standard deviations of the mean. As the NROY space is progressively reduced during successive waves, this threshold is lowered to 2 for the final wave. With the progressive reduction in emulator uncertainty, decreasing the implausibility cut-off is a standard step in the calibration procedure (Williamson et al., 2017; Hourdin et al., 2021).

# 3 New developments for polar snow in ORCHIDEE Explicit Snow

## 3.1 Snowpack initialization

Previously, ORCHIDEE was initialized without snow over ice sheet, relying on model accumulation through a five-year spin-up period (Charbit et al., 2024). Here, we introduce a new method for initializing the snowpack over the Antarctic and Greenland ice sheets. This procedure enables the formation of a deep snowpack in low accumulation areas, such as the interior of AIS where several years of simulation would be necessary to form a thick snowpack. Our initialization method expresses the initial density profile as a function of latitude and surface elevation only, making it easily adaptable to any snowpack model. This method is presented in details in Appendix A. Here, we summarize the main outlines. and summarized hereafter. To achieve this, we use the densification formulation proposed by Arthern et al. (2010) and further adapted by Ligtenberg et al. (2011):

$$\frac{d\rho}{dz} = C \cdot g \cdot \gamma(\ln(\bar{A})) \cdot \exp(\frac{E_g - E_c}{RT_s}) \cdot \rho(z) \cdot (\rho_{ice} - \rho(z)) \cdot dz \quad (2)$$

where $\bar{A}$ is the local mean annual accumulation (in kg m$^{-2}$ yr$^{-1}$), $T_s$ is the local mean annual surface temperature (in K), $g$ the gravitational acceleration (9.81 m s$^{-2}$), $\rho_{ice}$ the ice density (917 kg m$^{-3}$), $dz$ the snow layer thickness (in m), $R$ the gas constant (8.31 J K$^{-1}$ mol$^{-1}$), $E_g$ and $E_c$ the activation energies of grain growth (4.24 · 10$^4$ J mol$^{-1}$) and compaction (6 · 10$^4$ J mol$^{-1}$) respectively. $C$ is an unitless constant with different values above (0.03) and below (0.07) the critical density $\rho = 550$ kg m$^{-3}$, reflecting a faster densification rate near the surface.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

The $\gamma$ factor represents a correction to the accumulation term introduced by Ligtenberg et al. (2011). Their study revealed that the ratio of modeled to observed depths at which snow reaches critical densities of 550 kg m$^{-3}$ and 830 kg m$^{-3}$ is strongly correlated with mean annual accumulation across Antarctica. To account for this dependency, we use the same corrective term except that we applied a scaling factor of 0.9 for densities greater than 550 kg m$^{-3}$ to better match the observed density profile at a depth corresponding to 830 kg m$^{-3}$ in Antarctica:

$$ \gamma(x, y) = \begin{cases} \max(1.435 - 0.151 \cdot \ln(\bar{A}(x, y)), \gamma_{\min}) & \text{for } \rho \le 550 \text{ kg m}^{-3} \\ \max(0.9 \cdot [2.366 - 0.293 \cdot \ln(\bar{A}(x, y))], \gamma_{\min}) & \text{for } \rho \ge 550 \text{ kg m}^{-3} \end{cases} \quad (3) $$

where $\gamma_{\min}$ is the minimal correction equal to 0.25 as in Ligtenberg et al. (2011).

After integration of Eq. (2) (Appendix A), the explicit formulation of the density profile is expressed as a function of the mean surface temperature $\bar{T}_s$, the logarithm of the mean annual accumulation rate $\ln(\bar{A})$, and the surface density $\rho_s$. We estimate these key variables as functions of surface elevation and latitude, based on physical considerations and empirical relationships derived from observations and modeling outputs and detailed in Appendix A. Over Greenland, the parameterizations yield root mean square errors (RMSE) corresponding to 15% of the mean surface temperature (RMSE = 2.8 °C), 11% of the mean accumulation logarithm (0.60 $\ln(kg m^{-2} yr^{-1})$) and 11% of the mean surface density (36.4 kg m$^{-3}$). Over Antarctica, the corresponding errors represent 14% of the mean surface temperature (4.5 °C), 12.2% of the mean accumulation logarithm (0.56 $\ln(kg m^{-2} yr^{-1})$), and 8% of the mean surface density (29.4 kg m$^{-3}$). When evaluated over both ice sheets, the parameterized surface mean temperature correlates strongly with the observed 10 m snow temperature, with a bias of only 0.8 °C and a correlation coefficient of 0.93. (Fig. 2a). This initialization method can generate a polar snowpack of prescribed depth over the AIS and GrIS. As ablation is negligible in Antarctica, we consider a spatially uniform snowpack thickness set to 10 m. Conversely, in Greenland, the snow depth varies spatially from the edges to the interior of the ice sheet where it is also fixed at 10 m (see Appendix A). The estimated mean surface temperature of the snowpack is also used to initialize snow and ice temperatures at the beginning of the simulation, ensuring a consistent initial state across snow and ice layers.

The initialized snowpack density is compared to observations under dry snow conditions, in order to isolate the densification driven by the snow compaction (Fig. 2b, c). In ablation zones, where melt is significant, additional densification occurs from the refreezing of meltwater retained within the snowpack. This process is not accounted for in the initialization procedure. Overall, the initialized snowpack density agrees well with observations under dry snow conditions. The RMSE represents 7% of the mean snow density at 5 m depth (RMSE = 32.4 kg m$^{-3}$). At 10 m depth, the comparison with observations performs even better, with reduced RMSE (26.4 kg m$^{-3}$, 5.1% of the mean snow density at this depth) and bias (4.4 kg m$^{-3}$; 0.9%).
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 2.** (a) Mean surface temperature $\bar{T}_s$ obtained from the initialization procedure represented as a function of the observed 10 m snow temperature. (b, c) Initialized snow density at 5 m (b) and 10 m (c) depth versus observed densities at 5 m (b) and 10 m (c) across Greenland, East Antarctica, and West Antarctica. Statistical indicators (correlation coefficient $R^2$, RMSE, and bias) are reported in blue. Locations in Greenland, East Antarctica, and West Antarctica are marked in red, cyan, and green respectively.</figcaption></figure>

## 3.2 Surface snow density

Previous versions of the ORCHIDEE snow model (Wang et al., 2013; Charbit et al., 2024) represented surface snow density as a function of surface air temperature and wind speed, following Pahaut (1976). This formulation was developed for seasonal snow where fresh snow densities typically range from 50 to 150 kg m⁻³ (Valt et al., 2018). In contrast, over Greenland and Antarctica, surface snow densities measured in the uppermost centimetres of the snowpack are significantly higher, ranging from 200 to 400 kg m⁻³ (Fausto et al., 2018; Agosta et al., 2019). In this way, applied to GrIS and AIS this formulation leads to a significant underestimation of surface snow density when compared to observations (Section 4.1).

245

We propose a new parameterization of surface snow density, common to both ice sheets and depending only on wind speed. Here, we consider surface snow density as the mean density of the upper 20 cm. To represent high surface densities in polar regions, we choose to compensate the limited physical understanding of the processes involved, by introducing a surface density parameterization based only on wind speed. This formulation is applied to the modeled fresh snow density in ORCHIDEE, meaning that the freshly deposited snow directly represent surface snow density. Similar approaches have been employed in Arctic (Barrere et al., 2017) and Antarctic studies (Amory et al., 2021; Lenaerts et al., 2012; Keenan et al., 2021; Veldhuijsen et al., 2023). Moreover, using both Greenland and Antarctic surface snow observations, we find that the observed surface density correlates more strongly with wind speed than with air temperature, supporting our choice to use a wind parameterization for surface density (Fig. A3).

250

255

If surface density is expressed as a linear function of instantaneous wind speed, then the temporal average of surface density can be related to wind speed weighted by precipitation. In order to calibrate the formulation, we compare gridded-averaged surface snow density observations from Greenland and Antarctica with the precipitation weighted 10 m wind speed climatological average (1980-2019) simulated by MAR (Fig.3).

260
<figure><img src="image_15.png" /></figure>

In Antarctica, a linear relationship appears between wind speed and surface snow density. In contrast, Greenland shows a more heterogeneous pattern, with observations clustering around 320 kg m⁻³, in agreement with Fausto et al. (2018). To account for this heterogeneity we propose a new expression for surface snow density ρs, defined as a function of the wind speed at 10 m,

$$ \rho_s(\text{ws}_{10}) = \rho_{\text{s MIN}} + (\rho_{\text{s MAX}} - \rho_{\text{s MIN}}) \cdot (1 + \tanh(\frac{\text{ws}_{10} - \text{ws}_{10 \text{CUT}}}{\text{ws}_{10 \text{INT}}})) / 2 \quad (4) $$

with ρs MIN = 260 kg m⁻³ the minimum surface snow density reached at low wind speeds (<5 m s⁻¹), ρs MAX = 450 kg m⁻³ the maximum surface snow density reached at high wind speeds (>15 m s⁻¹), ws 10 CUT = 8.5 m s⁻¹ the wind speed at the inflection point of the hyperbolic tangent, ws s INT = 4.0 m s⁻¹ the width of the interval over which the hyperbolic transition from minimum to maximum density takes place. For Greenland, we imposed a maximum density of 400 kg m⁻³ to ensure a realistic range of surface densities in this region, as observations generally do not exceed this value (Fausto et al., 2018). For Antarctica, the maximum density remains equal to ρs MAX.

<figure><img src="image_16.png" /><figcaption>**Figure 3.** Observed mean surface snow densities in the upper 20 cm of the snowpack represented as a function of the climatological average (1980-2019) of the 10 m wind speed simulated by MAR and weighted by the amount of snowfall during the same period. Observations are presented in blue for Greenland and orange for Antarctica. The red line represents the wind speed parameterization of surface snow density presented in this study.</figcaption></figure>

Applied to both ice sheets, the formulation yields a RMSE of 40.3 kg m⁻³ (11.3% compared to the averaged surface density) and a bias of 1.9 kg m⁻³ (i.e. 0.5%), indicating a satisfactory agreement with observations.
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

## 3.3 Automatic calibration of snow densification formulations

### 3.3.1 Snow densification in ORCHIDEE

The parameterization of snow density in ORCHIDEE (Anderson, 1976) accounts for two processes: the destructive metamorphism and the compaction by overburden pressure. It represents dry snow densification and is expressed as follows:

$$ \frac{1}{\rho} \cdot \frac{\partial \rho}{\partial t} = \psi(T, \rho) + \frac{\sigma}{\eta(T, \rho)} \qquad (5) $$

where $\rho$ and $T$ are the snow density (in kg m⁻³) and temperature (in K) respectively, $\psi$ is the destructive metamorphism term (in s⁻¹), $\sigma$ represents the vertical stress of the overlying snow (in Pa) and $\eta$ is the snow viscosity (in kg s⁻¹ m⁻²).

Destructive metamorphism occurs primarily under isothermal or small temperature gradient conditions within the snow-pack, and influences directly the bonding between snow crystals. It refers to the movement of water molecules across the snow crystals through phase changes such as sublimation and condensation (van Kampenhout et al., 2017). The destructive metamorphism term is given by the following expression:

$$ \psi = c_{\psi} \cdot \exp(-c_{T} \cdot (T_{f} - T) - c_{\rho} \cdot \max(0, \rho - \rho_{\text{lim}})) \qquad (6) $$

where $T_f$ is the fusion temperature (273.15 K), $c_\psi = 2.8 \cdot 10^{-6} s^{-1}$, $c_T = 4.2 \cdot 10^{-2} K^{-1}$, $c_\rho = 460 \text{ m}^{-3} \text{ kg}^{-1}$ and $\rho_{\text{lim}} = 150 \text{ kg m}^{-3}$. This term is only active for densities lower than $\rho_{\text{lim}}$, and decreases exponentially below this threshold. In Greenland and Antarctica, surface snow densities exceed $\rho_{\text{lim}}$, implying that this term can be neglected in the rest of this study.

The second term of the right hand side of Eq. (5) represents the compaction by overburden pressure. It is driven by the vertical stress of the overlying snow, which induces sintering, reduces porosity, and leads to continuous deformation through mechanical creep. This process is the dominant contribution to dry snow densification. Anderson (1976) expressed the densification rate due to the overburden pressure as the second term of Eq. (5). The snow viscosity $\eta$, is expressed as a function of density and temperature:

$$ \eta = \eta_0 \cdot \exp[a_{\eta}(T_f - T) + b_{\eta}\rho] \qquad (7) $$

where $\eta_0 = 3.7 \cdot 10^7 \text{ kg s}^{-1} \text{ m}^{-2}$, $a_\eta = 8.1 \cdot 10^{-2} \text{ K}^{-1}$ and $b_\eta = 1.8 \cdot 10^{-2} \text{ m}^3 \text{ kg}^{-1}$.

The values of the snow viscosity parameters provided by Anderson (1976) are not adapted to polar snow conditions, leading to excessive densification compared to observed density profiles (van Kampenhout et al., 2017) (Fig. 1c,f).

### 3.3.2 Metrics for densification and uncertainty quantification

In this study, we use the HM method, described in Section 2.4, to calibrate three parameters of the snow viscosity defined by Anderson (1976): $\eta_0$, $a_\eta$ and $b_\eta$. We calibrate ORCHIDEE using 1D simulations over two dry snow locations, one in Greenland
<figure><img src="image_19.png" /></figure>

(Summit, 5 density profiles) and one in Antarctica (Dome C, 2 density profiles), presented in Section 2.1 (Fig. 1c, f). For each density profile, we apply a third degree polynomial interpolation (called 'fit' hereafter) to smooth the observations. In this way, in the 1D simulations, the surface snow density ($\rho_s$) is prescribed using the average of the fitted surface snow density (Table 2). In addition, to characterize the dry snow densification profiles, we define two metrics for each location: density at 5 m ($\rho_{5m}$) and at 10 m ($\rho_{10m}$). Each metric is obtained from the fit of the observed density profiles, and is associated with an uncertainty estimate that accounts for both observational and model-related uncertainty.

Observational uncertainties are estimated by combining two sources of uncertainty. We distinguish *within-profile* uncertainty, which captures the vertical variability for a given profile around the considered depth, and *between-profiles* uncertainty, which reflects variability across different profiles at the same site, related to either spatial or temporal variability. Here, we present the observational uncertainties for the metrics ($\rho_{5m}$ and $\rho_{10m}$) and the observed surface density ($\rho_s$). First, *within-profile* variability is quantified as the Root Mean Square (RMS) of the mean deviation (taken as $\frac{\sigma}{\sqrt{n}}$) across a vertical window around the target depth (4 to 6 m for $\rho_{5m}$, 9 to 11 m for $\rho_{10m}$, 0 to 1 m for $\rho_s$), e.g. for $\rho_{5m}$:

$$ \sigma_{within-profile 5m} = \sqrt{\frac{1}{l} \sum_{i=1}^{l} \left( \frac{\sigma(\rho_{obs i 4m<z<6m})}{\sqrt{n_i}} \right)^2} \quad (8) $$

where $l$ is the number of observed density profiles available, $n_i$ the number of observations of a given profile $i$ within the vertical window around the target depth, $\rho_{smoothed i}$ and $\rho_{obs i}$ the observed density of the given profile $i$ within the same vertical window.

Second, *between-profiles* variability is estimated as the standard deviation of the density metrics across all profiles for each location, divided by the square root of the number of profiles, e.g. for $\rho_{5m}$:

$$ \sigma_{between-profiles 5m} = \frac{\sigma(\rho_{smoothed 5m 1,...,l})}{\sqrt{l}} \quad (9) $$

where $\rho_{smoothed 5m 1,...,l}$ is the vector of the $l$ density metrics (here $\rho_{5m}$) resulting from the polynomial interpolation of the corresponding observed density profile.

These two components ($\sigma_{within-profile}$ and $\sigma_{between-profiles}$) are then combined using the root sum square to consider independently both uncertainties, yielding to the observational uncertainty for each metric ($\sigma_{obs 5m}$, $\sigma_{obs 10m}$) and for surface density ($\sigma_{\rho_s obs}$). In this way, at 5 m, mean densities from fitted profiles are equal to 409.1 and 367.3 kg m⁻³ at Summit and Dome C, associated with observational uncertainties equal to ±9.3 and ±5.7 kg m⁻³. At 10 m, we obtain respectively 497.5 and 415.6 kg m⁻³, with observational uncertainties of ±6.8 and ±3.4 kg m⁻³. The observed surface densities are equal to 311.4 ±10.2 kg m⁻³ at Summit and 316.4 ±9.2 kg m⁻³ at Dome C (Table 2).

In addition to the observational uncertainty, we need to account for modeling uncertainties. We evaluate the sensitivity of the metrics $\rho_{5m}$ and $\rho_{10m}$ to biases in the modeled snow temperature, accumulation and surface snow density using the steady-state formulation of Ligtenberg et al. (2011) (Eq. (2), Section 3.1). These biases are evaluated by considering the 40-year integrated accumulation and the final year of snow temperature simulated by ORCHIDEE in comparison with the in situ observations
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

presented in Section 2.1. We obtain temperature biases of +2.0 °C at Dome C and +2.6 °C at Summit, and accumulation biases of 12.4 kg m⁻² yr⁻¹ at Dome C and -54 kg m⁻² yr⁻¹ at Summit. We compute new steady state density profiles after adding the absolute value of these biases to the observed temperature and accumulation, using both positive and negative perturbations. The difference between the two perturbed steady state profiles is assumed to represent the model sensitivity to temperature and accumulation biases. The uncertainty is then estimated for ρ₅m and ρ₁₀m, as half this difference (δₐccu model and δT model). The same method is applied for ρs by injecting observed surface snow density uncertainty σobsρs, defined previously, in the steady state densification equation for each location separately.

Then, the final model uncertainty for each metric is obtained by Root Sum Square aggregation of the uncertainties presented previously, e.g. for ρ₅m :

$$ \sigma_{\rho_{5m}} = \sqrt{\sigma_{obs}^2 \rho_{5m} + \delta_{accu model 5m}^2 + \delta_{T model 5m}^2 + \delta_{\rho_s model 5m}^2} \quad (10) $$

This yields total uncertainties of ±15.3 kg m⁻³ and ±19.2 kg m⁻³ for ρ₅m and ρ₁₀m at Summit, and ±15.0 kg m⁻³ and ±14.0 kg m⁻³, respectively, at Dome C (Table 2).

<table><caption>**Table 2.** Observed metrics and associated uncertainties for surface density (ρs), densities at 5 m (ρ₅m) and 10 m (ρ₁₀m) for Summit and Dome C. Observational (σobs) and model uncertainties (δT model, δaccu model, δρs model) are reported, as well as the total uncertainties (σtot).</caption><thead><tr><th rowspan="2">Metric (kg m<sup>-3</sup>)</th><th colspan="3">Summit</th><th colspan="3">Dome C</th></tr><tr><th>ρ<sub>s</sub></th><th>ρ<sub>5m</sub></th><th>ρ<sub>10m</sub></th><th>ρ<sub>s</sub></th><th>ρ<sub>5m</sub></th><th>ρ<sub>10m</sub></th></tr></thead><tbody><tr><td></td><td>311.4</td><td>409.1</td><td>497.5</td><td>316.4</td><td>367.3</td><td>415.6</td></tr><tr><td>σ<sub>within-profile obs</sub></td><td>9.1</td><td>6.5</td><td>5.6</td><td>8.9</td><td>5.7</td><td>3.4</td></tr><tr><td>σ<sub>between-profile obs</sub></td><td>4.7</td><td>6.6</td><td>3.8</td><td>0.3</td><td>0.1</td><td>0.3</td></tr><tr><td>σ<sub>obs</sub></td><td>10.2</td><td>9.3</td><td>6.8</td><td>8.9</td><td>5.7</td><td>3.4</td></tr><tr><td>δ<sub>T model</sub></td><td></td><td>6.2</td><td>12.7</td><td></td><td>3.5</td><td>7.2</td></tr><tr><td>δ<sub>accu model</sub></td><td></td><td>3.8</td><td>7.8</td><td></td><td>2.9</td><td>6</td></tr><tr><td>δ<sub>ρs model</sub></td><td></td><td>9.8</td><td>10.1</td><td></td><td>9.5</td><td>9.8</td></tr><tr><td>σ<sub>tot</sub></td><td></td><td>15.3</td><td>19.2</td><td></td><td>15.0</td><td>14.0</td></tr></tbody></table>

### 3.3.3 Snow viscosity calibration with 1D simulations

The HM method is applied to calibrate the viscosity parameters of Eq. (5), with the objective of identifying parameter sets consistent with snow densification at both Summit and Dome C. Since we calibrate three viscosity parameters, the sample size is equal to 3 × 10 sets of parameters, sampled at each wave of the HM algorithm, following Loeppky et al. (2009) (see also Section 2.4). For each parameter set and each location, we run 1D ORCHIDEE simulations (1 at Summit and 1 at Dome C) and compute the corresponding metrics. It results in 60 simulated ρ₅m and ρ₁₀m values for each wave (3 parameters x 10 simulations x 2 locations), used to construct the Gaussian process emulators. In total, we perform 5 waves, to progressively reduce the parameter space and obtain parameter values providing consistent densification metrics ρ₅m and ρ₁₀m.
<figure><img src="image_22.png" /></figure>

During the first wave, the simulated densities at 5 m and 10 m depths exhibit a wide spread, ranging from 300 to 900 kg m⁻³, highlighting the strong sensitivity of modeled densification to the viscosity parameters (Fig. 4a, b). After five waves, the parameter space is reduced by 99.6%, and the resulting viscosity parameters are able to capture observed densification rates. At wave 5, the simulated ρ₅m and ρ₁₀m remain within the bounds of the metrics uncertainties at both Summit and Dome C.

<figure><img src="image_23.png" /><figcaption>**Figure 4.** (a, b) Evolution of density profiles simulated by ORCHIDEE in 2019, after 40 years of simulation for Summit (a) and Dome C (b) along the first and fifth waves of the HM method. Waves 1 and 5 are respectively represented in pink and blue. (c, d, e) Density metrics at 5 m and 10 m are shown with dots. Reduction of the parameter ranges of η₀ (c), a_η (d) and b_η (e) between the first and the fifth wave, compared to the RMSE accounting for the metrics ρ₅m and ρ₁₀m at each location.</figcaption></figure>
<figure><img src="image_24.png" /></figure>

The HM method also leads to a substantial reduction of the viscosity parameter ranges between the first and fifth wave (Fig. 4c, e; Table 3). The sensitivity of snow viscosity to temperature, governed by $a_\eta$, is increased from 0.081 to values converging between 0.138 and 0.17 K$^{-1}$, highlighting the stronger control of temperature on polar densification compared to the former formulation. This suggests that snow viscosity was underestimated in the original parameterization (Anderson, 1976), due to an insufficient influence of the temperature.

In contrast, the sensitivity of the viscosity to density, controlled by parameter $b_\eta$, remains close to the values proposed by Anderson (1976), with a final range between 0.018 and 0.02 m$^{-3}$ kg$^{-1}$ (0.018 m$^{-3}$ kg$^{-1}$initially). By comparison, the parameter $\eta_0$, which sets the reference viscosity, exhibits limited influence on the metrics, as its variability does not significantly affect model performance within the explored range. This suggests that either the selected metrics are not sensitive enough to constrain $\eta_0$ or that polar snow densification is primarily driven by the temperature and density sensitivities ($a_\eta$ and $b_\eta$), with $\eta_0$ playing only a minor role. The analysis of the final parameter space reveals a strong correlation between both $a_\eta$ and $b_\eta$ (Fig. A4). These parameter refinements are reflected in the reduction of the aggregated RMSE, i.e. the RMSE computed combining the four density metrics (at 5 m and 10 m depths over Summit and Dome C):

$$ \sqrt{\frac{1}{n} \sum_{i=0}^{n} [(\rho_{i5m} - \rho_{i5m \text{ obs}})^2 + (\rho_{i10m} - \rho_{i10m \text{ obs}})^2]} \quad (11) $$

where n is the number of locations studied (Summit and Dome C). The aggregated RMSE decreases substantially from an initial range of 50 to 500 kg m$^{-3}$ to a reduced range between 6.85 and 9.87 kg m$^{-3}$ after five calibration waves (Fig. 4c, d, e). We select the parameter set that minimizes the aggregated RMSE (6.85 kg m$^{-3}$) across both locations (Table 3).

<table><caption>**Table 3.** Free parameters of the Anderson (1976) viscosity equation used in the calibration. “Anderson value” corresponds to the default value of each parameter, “Prior range” to the initial parameter bounds, and “Calibrated range” and “Calibration result” to the ranges and values retained after the HM method.</caption><thead><tr><th>Parameter</th><th>Controls</th><th>Anderson value</th><th>Prior Range</th><th>Calibrated Range</th><th>Calibration result</th></tr></thead><tbody><tr><td>η<sub>0</sub> [kg s<sup>-1</sup> m<sup>-2</sup>]</td><td>Reference viscosity</td><td>3.7 &middot; 10<sup>7</sup></td><td>[3 &middot; 10<sup>7</sup>, 7.7 &middot; 10<sup>7</sup>]</td><td>[3.06 &middot; 10<sup>7</sup>, 6.91 &middot; 10<sup>7</sup>]</td><td>5.0216 &middot; 10<sup>7</sup></td></tr><tr><td>a<sub>η</sub> [K<sup>-1</sup>]</td><td>Viscosity dependency to temperature</td><td>0.081</td><td>[0.001, 0.25]</td><td>[0.138, 0.17]</td><td>0.15297</td></tr><tr><td>b<sub>η</sub> [m<sup>-3</sup> kg<sup>-1</sup>]</td><td>Viscosity dependency to ρ</td><td>0.018</td><td>[0.001, 0.03]</td><td>[0.018, 0.02]</td><td>0.01952</td></tr></tbody></table>

## 4 Application to the Greenland and Antarctic ice sheets

The developments presented in the previous sections have led to noticeable improvements in the simulated density profiles for Summit and Dome C compared to available observations. In this section, we now perform offline ORCHIDEE simulations (1980-2019) for the entire GrIS and AIS to estimate their response to these new developments.
<figure><img src="image_25.png" /></figure>

## 4.1 Surface snow density

We compare the new formulation of the surface snow density computed by ORCHIDEE with observations of surface snow densities (Fig. 5). This comparison is made for the former formulation (Pahaut, 1976) and the new one proposed in this study (Eq. 4). Surface snow densities are computed as the average density of the top 20 cm of the snowpack over a 40 year period (1980-2019).

As previously mentioned, the Pahaut (1976) parameterization strongly underestimates surface snow density for both AIS and GrIS, with a negative bias of -188.3 kg m⁻³, representing 53 % of the average observed ρs. The new formulation shows an overall good representation of surface density with a bias of only 5.8 kg m⁻³. Regional analysis reveals a positive bias over Greenland (20.6 kg m⁻³, i.e. 6.2 % with reference to the observed average ρs) and a negative bias over Antarctica (-23.8 kg m⁻³, 6.3 %). However, surface snow densities can be affected by melt-refreeze processes, particularly over Greenland. These processes increase snow density through the accumulation of refrozen meltwater. Furthermore, this additional densification is represented within the snow scheme, but is not explicitly included in the surface snow density parameterization. This can constitute a potential bias for our evaluation.

We account for this effect by excluding simulated surface snow densities where the 40 year average refreezing rate in the top 20 cm of the snowpack exceeds 0.1 kg m⁻² d⁻¹. These locations are represented by transparent dots in Figure 5a.

<figure><img src="image_26.png" /><figcaption>**Figure 5.** (a) Comparison of 40 year averaged (1980-2019) modeled surface snow density (computed as the average density of the top 20 cm of the snowpack) with observed surface snow density from Greenland (open diamond) and Antarctica (solid circle). Model results are shown for the previous ORCHIDEE formulation (OR OLD, purple; Pahaut, 1976) and the new formulation (OR NEW, blue). (b, c) Spatial comparison of 40-year averaged modeled surface snow density with observations (in dots) in Antarctica (b) and Greenland (c).</figcaption></figure>

Excluding the Greenland locations affected by refreeze processes improves the computed surface density, when considering both ice sheets, with a negative bias of -3.5 kg m⁻³ and a reduced RMSE equal to 34.8 kg m⁻³ (1 % and 9.8 % compared to the
<figure><img src="image_27.png" /></figure>

observed averaged surface snow density, respectively). This highlights the need for future studies to investigate the ORCHIDEE representation of refrozen areas in Greenland.

Moreover, the new parameterization generally reproduces the large-scale spatial gradients in surface snow density over both ice sheets (Figs. 5b,c). It captures the higher densities along the wind exposed coastal margins, especially over Antarctica.

In addition, we conduct an inter-model comparison by implementing the surface snow parameterizations from both MAR (Kittel et al., 2021) and IMAU-Firn densification model (IMAU-FDM; Lenaerts et al., 2012; Veldhuijsen et al., 2023) within ORCHIDEE. The MAR formulation of surface snow density is based only on 10 m wind speed and was developed for Antarctica. The IMAU-FDM formulation relies on the 10 m wind speed and near surface temperature ($T_{2m}$) and was also calibrated for Antarctic conditions. We test these formulations with offline ORCHIDEE simulations over Greenland and Antarctica for the 1980–2019 period. Model outputs are evaluated considering the top 20 cm for surface snow densities (Fig. 6). Both the IMAU-FDM and MAR parameterizations perform well over Antarctica, with biases of only -1.3 and 1.5 kg m⁻³ respectively, representing differences of only 0.3 % and 0.4 % with the averaged observed surface density (Table 4). However, their performance degrades over Greenland, where biases increase to 60.8 kg m⁻³ (18.4 %; IMAU-FDM) and 65.3 kg m⁻³ (19.8 %; MAR). This behavior highlights the challenge of representing both ice sheets with a single parameterization, given their distinct surface conditions. Our proposed formulation reproduces reasonably surface density in both regions, with biases of only 6.2 % and 7.2 % (20.6 kg m⁻³ and -23.8 kg m⁻³) for Greenland and Antarctica respectively, and a relatively small RMSE (47.3 kg m⁻³ in Greenland; 44.1 kg m⁻³, in Antarctica). Combined over both ice sheets, the proposed formulation gives a reduced bias (5.8 kg m⁻³, 1.6 %) and RMSE (44.1 kg m⁻³, 12.4 %).

<table><caption>**Table 4.** Surface snow density statistical indicators for MAR, IMAU and ORCHIDEE surface snow formulations</caption><thead><tr><th rowspan="2"></th><th colspan="3">Greenland</th><th colspan="3">Antarctica</th><th colspan="3">Both IS</th></tr><tr><th>RMSE</th><th>BIAS</th><th>R<sup>2</sup></th><th>RMSE</th><th>BIAS</th><th>R<sup>2</sup></th><th>RMSE</th><th>BIAS</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><th>MAR</th><td>77.7</td><td>65.3</td><td>0.11</td><td>30.5</td><td>1.5</td><td>0.43</td><td>65.9</td><td>44.2</td><td>0.05</td></tr><tr><th>IMAU</th><td>75.5</td><td>60.8</td><td>0.11</td><td>29.8</td><td>-1.3</td><td>0.55</td><td>64.1</td><td>40.2</td><td>0.1</td></tr><tr><th>OR NEW</th><td>47.3</td><td>20.6</td><td>0.09</td><td>36.9</td><td>-23.8</td><td>0.56</td><td>44.1</td><td>5.8</td><td>0.18</td></tr></tbody></table>
<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 6.** 40y-average surface snow density compared to observations, modeled by (a) the new wind parameterization in ORCHIDEE, (b) the wind parameterization developed for Antarctica in MAR (Kittel et al., 2021) and (c) the wind and temperature parameterization developed for Antarctica in IMAU (Lenaerts et al., 2012). Locations are specified with open diamond and solid circle, for Greenland and Antarctica respectively. Shaded points are locations where the 40y-averaged refreeze in the top 20 cm of the snowpack is superior to 0.1 mm d⁻¹</figcaption></figure>

## 4.2 Polar snow densification in ORCHIDEE

To assess whether the densification improvements obtained in 1D simulations extend to the ice sheet scale, we now compare modeled and observed snow densities at 5 m and 10 m in 2D configurations for both ice sheets (Fig. 7).

The agreement between ORCHIDEE outputs and observations is generally good, with a negative bias of -2.3 kg m⁻³ (i.e. 0.5% of the ρ₅m average) and a RMSE of 42.8 kg m⁻³ (9.4% with reference to the average ρ₅m). A modest underestimation of density can be noted over Antarctica (Fig. 7a). A similar performance is obtained at 10 m with differences of 1.6% and 8.5% for the bias and the RMSE respectively. These statistics demonstrate that the model reproduces the vertical structure of observed density profiles with reasonable accuracy. To further illustrate the model performance in reproducing the ρ₅m and ρ₁₀m observed metrics, we perform an evaluation over all available locations (see Figs. S1, S2 and S4). Here, we selected ten dry snow observation sites (five in Greenland and five in Antarctica), chosen to cover the different regions of both ice sheets. For each location, the previous ORCHIDEE configuration is compared with the updated version of the model, including the new developments (Fig. 7b). Density profiles correspond to the annual average for the year 2019, computed from monthly outputs. Two systematic biases appear consistently across all sites when we consider the previous representation of density in ORCHIDEE (OR OLD). First, the model underestimates surface snow density by about 150 kg m⁻³. Second, it overestimates the compaction, leading to an increased density within the snowpack of around 100 to 200 kg m⁻³ compared to the observations. In the updated version of ORCHIDEE (OR NEW), these biases are considerably reduced and the vertical structure of snow density is significantly improved, despite a slight underestimation in some specific Antarctic locations such as Mount Johns (Fig. 7b).
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 7.** Comparison of (a) densities at 5 m and 10 m depth between observations and the ORCHIDEE model with the calibrated set of snow viscosity parameters across both GrIS and AIS. (b) Ten density profiles representative of dry snow locations across both ice sheets. Dark lines represent observed density profiles, orange dashed lines indicate the initialization profile, the purple curves (OR OLD) show the density profile before developments and calibration and the blue ones (OR NEW) represent the new density profile.</figcaption></figure>

## 4.3 Effect of initialization and model developments on snow density

To investigate the effect of the initialization procedure as well as the developments proposed in this study, we choose to analyze the monthly evolution of the modeled snow density at Summit under three different ORCHIDEE configurations (Table 5, Figure 8): (i) without initialization meaning that we start the simulation from a snow-free state over the ice sheet (OR OLD),

435
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

(ii) with standard initialization (OR OLD+INIT), and (iii) with initialization combined with updated model physics, including calibrated viscosity parameters for polar regions and new surface density parameterization (OR NEW).

<table><caption>**Table 5.** Overview of the model configurations and developments included in each simulation. OR OLD corresponds to the original model configuration, OR OLD+INIT includes the snowpack initialization procedure, and OR NEW takes into account the different developments presented in this study.</caption><thead><tr><th>Category</th><th>Developments</th><th>OR OLD</th><th>OR OLD+INIT</th><th>OR NEW</th></tr></thead><tbody><tr><td>Refreezing (Section 2.2)</td><td>Refreeze correction</td><td>×</td><td>×</td><td>×</td></tr><tr><td>Initialization (Section 3.1)</td><td>Snowpack initialization</td><td></td><td>×</td><td>×</td></tr><tr><td>Surface density (Section 3.2)</td><td>Pahaut (1976) parameterization<br>Wind-based parameterization</td><td>×</td><td>×</td><td>×</td></tr><tr><td>Snow compaction (Section 3.3)</td><td>Anderson (1976): initial viscosity parameters<br>Anderson (1976): calibrated viscosity parameters</td><td>×</td><td>×</td><td>×</td></tr></tbody></table>

In OR OLD, it takes approximately 30 years for the model to accumulate 10 meters of snow at Summit (Fig. 8a). The initialized simulation OR OLD+INIT begins with a 10 meter snowpack, with an initialized steady-state density consistent with observed density profiles at Summit (Fig. 8b,h). Rapid thinning of the snowpack is observed in the first few years due to the excessive densification rate that causes OR OLD to deviate quickly from the initial density profile.

We further analyze the convergence time of the model towards its equilibrium state by considering anomalies between monthly snow density and the average profile from the final simulation year (Fig. 8d,e,f). In OR OLD without initialization, it takes more than 35 years to converge towards the final density profile (Fig. 8d), while OR OLD+INIT converges slightly faster, within roughly 30 years (Fig. 8e). However, OR OLD+INIT clearly diverges from the initial profile during the first decade, as indicated by the pronounced negative anomaly during this period. By the end of the simulation, the resulting density profile is significantly denser by more than 150 kg m⁻³ compared to the initialized profile, particularly in the deeper layers.

Finally, when comparing monthly modeled density profiles with observed density profile at Summit (taken as the polynomial fit on the 05-2017 profile), both OR OLD and OR OLD+INIT produce too low densities near the surface (Fig. 8g,h), confirming known biases in surface densification presented in Section 4.1. Below the surface, the OR OLD densities are systematically too high compared to observations, regardless of the initialization (Fig. 8g,h). This confirms the excessive densification found by van Kampenhout et al. (2017) using the Anderson densification parameterization in polar regions.

In contrast, the OR NEW configuration exhibits a more realistic vertical structure, with densities ranging between 300 kg m⁻³ and 600 kg m⁻³ over the top 17 meters, closely matching the observed density profile (Fig. 8c,i). The bias remains below 30 kg m⁻³ through the snow column. Furthermore, OR NEW reaches equilibrium within approximately five years, substantially reducing the spin-up time (Fig. 8f). The anomaly with respect to the initial profile remains low during the entire simulation, confirming the consistency of the initialization with the long-term behavior of the model.
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /><figcaption>**Figure 8.** Monthly evolution of (a, b, c) the modeled snow density at Summit; (d, e, f) the density anomaly between modeled snow density and the final year averaged (2019) simulated profile, interpolated over the model snow depth at Summit for each month; and (g, h, i) the density anomaly between monthly simulated profiles and an observed density profile derived from a third-degree polynomial fit. The results are shown for three model configurations: (left) OR OLD, the previous version of ORCHIDEE, (middle) OR OLD+INIT with the initialization procedure activated, and (right) OR NEW, the final configuration with the initialization combined with the densification developments.</figcaption></figure>
<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

# 5 Impact of developments on temperature and melt-refreeze processes

## 5.1 Evaluation of the 10 m snow temperature

In the ORCHIDEE snow scheme, snow thermal conductivity and snow heat capacity are expressed as functions of density (Charbit et al., 2024). As a result, improvements in the representation of snow density should translate into a better representation of heat diffusion within the snowpack, and thus of the vertical temperature profile. In the following, the simulated temperature at 10 m ($T_{10m}$), yearly averaged over the year 2019, are compared to observations using both the previous configuration (OR OLD) and the updated version (OR NEW) including all developments. In the original configuration, $T_{10m}$ is persistently overestimated across both ice sheets, with a bias of 4.0 °C in Greenland and 3.1 °C in Antarctica (Fig. 9a,c). The updated version improves the performance over Greenland and Antarctica, reducing the bias to 1.5 °C and 1.4 °C respectively (Fig. 9b,d).

<figure><img src="image_39.png" /><figcaption>**Figure 9.** Comparison of observed and simulated snow temperature at 10 m depth from (a, c) the previous ORCHIDEE configuration (OR OLD) for (a) Greenland and (c) Antarctica; and from the new ORCHIDEE configuration (OR NEW) for (b) Greenland and (d) Antarctica.</figcaption></figure>
<figure><img src="image_40.png" /></figure>

## 5.2 Evaluation of melt-refreeze processes

To assess the influence of our developments on modeled melt-refreezing processes in Greenland, we focus our analysis on the southeastern coast of GrIS, a region characterized by high snow accumulation and intense ablation (Fig. 10a). This area is particularly relevant for investigating refreezing dynamics, with MAR modeled mean melt rates exceeding 4.5 m yr⁻¹ water equivalent (w.e.) over 2010–2019 and recurrent episodes of strong melting and subsequent refreezing (Fig. 10b,c). In Charbit et al. (2024), two sets of albedo parameters were tested, one from Raoult et al. (2023) (presented in Section 2.3) and one from a manual tuning approach, designed to obtain a better agreement of the simulated runoff between MAR and ORCHIDEE (OPT12L). As the parameters from Raoult et al. (2023) lead to an underestimation of melt-refreeze processes (Charbit et al., 2024), we perform the same simulations as OR NEW and OR OLD+INIT with the OPT12L albedo parameters, to verify that the influence of the developments is consistent with both albedo parameters. Overall, these sets of parameters lead to similar conclusions. Moreover, since our analysis focuses on melt-refreeze processes, we will use the OPT12L configuration in the following to facilitate the interpretation (Fig. 10). First we find that even with the OPT12L albedo parameters, ORCHIDEE underestimates both melt and refreezing fluxes in comparison to MAR (Fig. 10b,c). The mean seasonal melt modeled by MAR, averaged over the region presented in Fig. 10a, ranges approximately between 6.5 m yr⁻¹ w.e and almost 10 m yr⁻¹ w.e. between 2010 and 2019, while ORCHIDEE shows differences of about 2 to 3 m yr⁻¹ w.e.. Refreezing is also underestimated, with differences between OR OLD+INIT and MAR of about 0.5 to 1 m yr⁻¹ w.e. (Fig. 10c). Nevertheless, the updated version (OR NEW) simulates enhanced refreezing relative to the initial configuration (OR OLD+INIT).

To better understand these differences, we analyze the 2012 melt season, which featured exceptional melt events across GrIS (Fig. 10d,e). Although simulated melt rates are similar between OR OLD+INIT and OR NEW, refreezing diverges notably early in the melt season (May–June 2012). We further investigate the intense melt event of 3–14 May 2012 (highlighted in light blue in Fig. 10d,e) at a representative site (Fig. 10a). The modeled snow density, temperature, and refreezing concentration (computed as the refreezing content divided by the mass of each snow layers) are compared between nighttime (01:30–04:30 UTC) and daytime (16:30–19:30 UTC), corresponding approximately to local midnight and early afternoon, respectively (Fig. 10f-h).

While melt rates remain stable between OR OLD+INIT and OR NEW (from 4.8 m yr⁻¹ w.e. in OR OLD+INIT to 4.4 m yr⁻¹ w.e. in OR NEW in average during this melt event), the refreezing rate is twice in OR NEW compared to OR OLD+INIT (from 1.4 m yr⁻¹ w.e. to 3 m yr⁻¹ w.e. in OR NEW). Indeed, during daytime, the refreezing in OR NEW is greater and tends to occur deeper within the snow column than in OR OLD+INIT (Fig. 10h). In addition, OR NEW appears to have a colder internal snow temperature profile in comparison with OR OLD+INIT, during both day and nighttime (Fig. 10g). Between 2 and 4 m depth, OR NEW remains up to 2 °C colder than OR OLD+INIT. Because refreezing depends on the available thermal energy, the colder snow in the deeper layers of OR NEW promotes more efficient refreezing of percolated meltwater (Fig. 10h), particularly during daytime when meltwater is produced. Moreover, these snow temperature differences between OR NEW and OR OLD+INIT can be explained by the changes in snow densities between both simulations (Fig. 10f). Indeed, at depth, OR NEW exhibits lower snow densities than OR OLD+INIT. As snow thermal conductivity is parameterized
<figure><img src="image_41.png" /></figure>

as a function of density (5.1), this results in a reduced snow conductivity in OR NEW. Consequently, the snowpack can act as a more effective insulator in the depth of the snowpack, explaining the colder internal snow temperature profile obtained in OR NEW.

In conclusion, the new developments promote higher refreezing rates by modifying the snow density and viscosity profiles. In OR OLD+INIT, the snowpack remains almost too warm to refreeze additional meltwater, either during day or night. In contrast, the colder snowpack in OR NEW allows for greater refreezing potential, which is consistent with the need to increase the melt and refreeze content in line with MAR.

<figure><img src="image_42.png" /><figcaption>**Figure 10.** (a) Average melt from the MAR model over the southeast coast of the GrIS, where melt exceeds a threshold of 0.2 m yr⁻¹. (b–c) Mean seasonal melt and refreezing rates (m yr⁻¹), averaged over the region shown in (a). Results are presented for OR OLD+INIT (orange), OR NEW (blue), and MAR (dashed lines). The year 2012 is highlighted in light blue. (d–e) Seasonal three-hourly melt and refreezing rates during 2012, averaged over the same region; the second major melt event (3–14 May 2012) is shaded in blue. (f–h) Mean snow density (f), temperature (g), and refreezing concentration (h) during this second melt event, at the site indicated in (a). Nighttime (01:30–04:30 UTC) and daytime (16:30–19:30 UTC) averages are shown as dark and light colors, respectively. Refreezing concentration is computed as the refreezing content of each snow layer divided by its corresponding snow mass. The averaged melt in m yr⁻¹ during this second melt event at this location is indicated for both ORCHIDEE simulations.</figcaption></figure>

510

515
<figure><img src="image_43.png" /></figure>

## 5.3 Densification in the ablation areas

As shown previously, ORCHIDEE underestimates refreezing within the snowpack, linked to an underestimation of meltwater production. These biases have an impact on densification in the Greenland ablation zone, defined generally as the low elevation area of the ice sheet where annual snow melt exceeds accumulation. Indeed, the resulting meltwater significantly influences water percolation and thus refreezing within the snowpack, which in turn affects densification.

In the ablation zone, ORCHIDEE underestimates snow density at both 5 m and 10 m (Fig. 11a), with a larger bias at 5 m (-56.8 kg m⁻³) than at 10 m (-21.3 kg m⁻³). At 10 m, the influence of refreezing on densification diminishes, as the water does not percolate significantly at this depth, and the compaction becomes the dominant densification process. Since the compaction scheme was calibrated using dry snow conditions, where refreezing is negligible, the model performs more accurately at 10 m. These results likely indicate that accurate representation of density in ablation zones may depend on the ability of the model to simulate melt and refreeze processes.

Further qualitative analysis is done by comparing observed and modeled density profiles from ORCHIDEE in the ablation zone (Fig. 11b). The following profiles show some successful reproduction of densification trends in the ablation zone. Exception is made for the first profile, which represents a modeled snowpack without clear evidence of refreezing (with surface densities around 600 kg m⁻³), whereas density observations suggest that refreezing occurs, leading to a densification of the snowpack. Other available observed profiles can be found in Supplementary (Fig. S1, S2 and S3) and show similar results.

<figure><img src="image_44.png" /></figure>
<figure><img src="image_45.png" /></figure>

**Figure 11.** Observed density at 5 m ρ<sub>5m</sub> (a) and 10 m ρ<sub>10m</sub> (b) compared to modeled ones with the new viscosity parameters across the Greenland ablation locations. (c) Locations of the five profiles selected (d) Density profiles representatives of ablation locations across Greenland. Dark lines represent observed density profiles (OBS), orange dashed line the initialization profile (INIT), the purple curve the ORCHIDEE density profile before developments (OR OLD) and the blue curve with the new ORCHIDEE configuration (OR NEW). Locations of the five density profiles are shown in (c).

# 6 Discussion and conclusion

We have introduced a series of model developments in ORCHIDEE aimed at improving the representation of polar snow density over the Greenland and Antarctic ice sheets. The model performance has been evaluated through offline experiments and compared against in situ observations. We first proposed a new initialization procedure for the snowpack over polar ice sheets. This method allows to start the simulations with a prescribed snow depth and a vertical density structure consistent with field observations. This method induces a substantial reduction in spin-up duration, which is especially valuable for coupled atmosphere-surface simulations. By relying only on latitude and surface elevation to reconstruct the density profile, this approach is aimed to be generalizable to any polar snowpack model.

In addition, we introduced a new empirical parameterization of surface snow density for polar regions, valid for both ice sheets and based on the 10 m wind speed. It has been calibrated using the wind speed weighted by precipitation against surface density observations. This formulation is a pragmatic approach to represent the rapid surface densification processes that are not explicitly represented in the model. Enhancing wind densification is considered a temporary solution, pending a better understanding and representation of the physical processes controlling surface snow densification. As an example, Domine et al. (2016) highlighted that the upward transport of water vapor driven by strong temperature gradients in polar snowpacks can contribute to an increase in surface density, a process not captured by ORCHIDEE. To compensate for unresolved processes, similar wind-based strategies have also been adopted in other regional models, such as MAR (Kittel et al., 2021) and the IMAU-Firn densification model (IMAU-FDM; Lenaerts et al., 2012; Veldhuijsen et al., 2023). Moreover, the comparison with surface snow density observations shows that the proposed parameterization successfully corrects the large underestimation present over polar regions in the previous model version.

Finally, we calibrated the viscosity parameters from the Anderson (1976) densification equation (Eq. 5). Using the History Matching procedure, we achieved a significantly improved representation of snow compaction in dry-snow conditions. These new viscosity parameters enable to reduce the overestimated compaction previously identified in ORCHIDEE in various locations of the GrIS and AIS. They also contribute to a reduction of the snow temperature bias in the model. However, this improvement still does not fully capture densification in ablation areas, where refreezing of percolated liquid water drives additional densification, especially in the first meters of the snowpack. As discussed in Section 5.2, both meltwater production and refreezing remain underestimated in ORCHIDEE. Considering the influence of the liquid water retention capacity formulation (following the approach of Lafaysse et al. (2017)) may be an effective step towards a better representation of liquid water retention and refreezing processes. The current formulation of liquid water holding capacity increases for densities above 200 kg m⁻³ (Boone and Etchevers, 2001), meaning that a denser snow can retain more liquid water. However, alternative formulations (Lafaysse et al., 2017) simulate the opposite behavior, consistent with laboratory experiments (Coléou and
<figure><img src="image_46.png" /></figure>

Lesaffre, 1998). Exploring the sensitivity of the model to these different parameterizations may help assess their impact on the simulation of refreezing processes.

Another perspective to improve the simulation of meltwater production, and therefore refreezing, is to focus on the representation of the surface energy balance, with particular attention to the snow albedo parameterization. This perspective will be explored using coupled experiments in order to capture surface/atmosphere feedbacks, particularly important for the albedo modeling. Special emphasis will be placed on refining the current simplified albedo formulation (Chalita and Le Treut, 1994), which represents albedo as a time-evolving function of snow age. Yet, the snow microstructure has a significant influence on snow spectral albedo (Skiles et al., 2018). This underlines the need for an explicit representation of the specific surface area in order to have a more physically based albedo scheme. Pending this development, an alternative approach consists of introducing different albedo decay rates for dry and wet snow in the current albedo parameterization, inspired by Helsen et al. (2017). This can be a step towards a better representation of the albedo decay rate during the melt season. Nevertheless, the developments presented in this paper improve significantly the representation of snowpack internal processes and are a key milestone for a more accurate assessment of the surface mass balance of polar ice sheets in the IPSL climate model.

## Appendix A: Snowpack Initialization

We propose an explicit formulation for the density $\rho$ after integration of the Ligtenberg et al. (2011) equation along the vertical axis, expressed as a differenced function based on the critical depth $z_{550}$ where $\rho$ equals 550 kg m$^{-3}$:

$$ \rho = \frac{\exp(E_0(\bar{T}_s)z)}{\frac{\rho_i}{\rho_s} - 1 + \exp(E_0(\bar{T}_s)z)} \rho_i \quad \text{if } z < z_{550} \qquad (A1) $$

$$ \rho = \frac{\exp(E_1(\bar{T}_s)(z - z_{550}))}{\frac{\rho_i}{z_{550}} - 1 + \exp(E_1(\bar{T}_s)(z - z_{550}))} \rho_i \quad \text{if } z > z_{550} \qquad (A2) $$

Here, $\rho_s$ refers to the surface density in kg m$^{-3}$, $E_0(\bar{T}_s) = C_0 g \exp(E_1(\bar{T}_s)) \gamma$ and $E_1(\bar{T}_s) = C_1 g \rho_i \exp(E_1(\bar{T}_s)) \gamma$, where $\gamma$ is the correction term defined previously in Section 3.1, and $C_0$ and $C_1$ are the constants introduced in Equation (1).

To derive a density parameterization that depends only on latitude and elevation, we aim to parameterize three key variables required in formulations (A1) and (A2) : the average surface temperature $\bar{T}_s$, the logarithm of the average annual accumulation rate $\ln(\bar{A})$, and the surface density $\rho_s$. These parameterizations are based on physical considerations and empirical relationships between these variables and surface elevation and latitude.

First, we introduce a parameterization for the mean surface temperature as a function of latitude and elevation, following the elevational lapse rates identified by Feulner et al. (2013) for Antarctica and by Hanna et al. (2005) for Greenland. Next, we define a parameterization for the logarithm of the mean accumulation rate $\ln(\bar{A})$ as a function of the mean annual surface
<figure><img src="image_47.png" /></figure>

temperature, based on the Clausius-Clapeyron relationship. Finally, we establish a parameterization for the surface density assuming a linear relationship between the surface density and the mean annual surface temperature. Each parameterization is defined separately for Greenland and Antarctica.

## Parameterization of mean surface temperature $\overline{T}_s$

To parameterize the mean surface temperature, we base our approach on the theoretical temperature gradients with respect to altitude and latitude for Antarctica (AIS) and Greenland (GrIS).

For the AIS, we adopted the dry-adiabatic lapse rate with respect to altitude, proposed by Feulner et al. (2013) and equal to $9.8^\circ\text{C km}^{-1}$. In our approach, we distinguished the temperature gradient with elevation and latitude based on a critical altitude of 500 meters above sea level ($h_{500}$). This threshold distinguishes two distinct regimes observed in the tendencies of the 10-meter snow temperature in the AIS (Fig. A2): a first tendency for elevation lower than 500 m a.s.l. taking into account ice shelves and the second regime associated to the rest of the ice sheet. For the GrIS, we followed the theoretical lapse rate with respect to elevation proposed by Hanna et al. (2005) equal to $5.9^\circ\text{C km}^{-1}$ below 1000 m a.s.l. ($h_{1000}$) and $8.2^\circ\text{C km}^{-1}$ above this critical altitude. The dependency of the mean surface temperature to latitude is based on an empirical relationship in order to represent the decrease of temperature towards the pole. To determine this temperature gradient, we applied the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization algorithm to minimize the root mean square error between the parameterized mean surface temperature and the observed firn temperature at a depth of 10 meters. We obtain the following temperature gradients with respect to latitude, resulting from the optimization and with respect to elevation based on the considerations of Feulner et al. (2013) and Hanna et al. (2005). For AIS:

$$ \left\{ \begin{array}{ll} \frac{\Delta T}{\Delta \text{lat low alt}} = 1.14^\circ\text{C lat}^{-1}, & \frac{\Delta T}{\Delta z \text{ low alt}} = 6.5^\circ\text{C km}^{-1} & \text{if } z < h_{500} \text{ m a.s.l,} \\ \\ \frac{\Delta T}{\Delta \text{lat high alt}} = 0.73^\circ\text{C lat}^{-1}, & \frac{\Delta T}{\Delta z \text{ high alt}} = 9.8^\circ\text{C km}^{-1} & \text{if } z \ge h_{500} \text{ m a.s.l.} \end{array} \right. \quad (A3) $$

For GrIS:

$$ \left\{ \begin{array}{ll} \frac{\Delta T}{\Delta \text{lat low alt}} = -0.66^\circ\text{C lat}^{-1}, & \frac{\Delta T}{\Delta z \text{ low alt}} = 5.9^\circ\text{C km}^{-1} & \text{if } z < h_{1000} \text{ m a.s.l,} \\ \\ \frac{\Delta T}{\Delta \text{lat high alt}} = -0.92^\circ\text{C lat}^{-1}, & \frac{\Delta T}{\Delta z \text{ high alt}} = 8.2^\circ\text{C km}^{-1} & \text{if } z \ge h_{1000} \text{ m a.s.l.} \end{array} \right. \quad (A4) $$

We express the parameterization of mean annual surface temperature as follows, depending on the temperature gradient with respect to latitude and elevation:
<figure><img src="image_48.png" /></figure>

$$ \frac{\Delta T}{\Delta \text{lat}}(z) = \left( \frac{\frac{\Delta T}{\Delta \text{lat high alt}} - \frac{\Delta T}{\Delta \text{lat low alt}}}{z_{\text{high alt}} - z_{\text{low alt}}} \right) \cdot (z - z_{\text{low alt}}) + \frac{\Delta T}{\Delta \text{lat low alt}} \quad (A5) $$

$$ \frac{\Delta T}{\Delta z}(z) = \left( \frac{\frac{\Delta T}{\Delta z \text{ high alt}} - \frac{\Delta T}{\Delta z \text{ low alt}}}{z_{\text{high alt}} - z_{\text{low alt}}} \right) \cdot (z - z_{\text{low alt}}) + \frac{\Delta T}{\Delta z \text{ low alt}} \quad (A6) $$

$$ T_s(\text{lat}, z) = \min \left( \frac{\Delta T}{\Delta \text{lat}} \cdot (\text{lat} - \text{lat}_{\text{ref}}) - \frac{\Delta T}{\Delta z} \cdot \frac{(z - z_{\text{ref}})}{1000} + T_{s\_\text{ref}}, 0 \right) \quad (A7) $$

The parameters used in the surface temperature parameterization for Antarctica and Greenland can be found in the following table.

<table><caption>**Table A1.** Parameters used in the surface temperature parameterization</caption><thead><tr><th>Parameter</th><th>Greenland</th><th>Antarctica</th></tr></thead><tbody><tr><td>T<sub>s_ref</sub> [°C]</td><td>1.11</td><td>0</td></tr><tr><td>lat<sub>ref</sub></td><td>60</td><td>-57</td></tr><tr><td>z<sub>ref</sub> [m]</td><td>1000</td><td>0</td></tr><tr><td>z<sub>low alt</sub> [m]</td><td>0</td><td>0</td></tr><tr><td>z<sub>high alt</sub> [m]</td><td>1000</td><td>500</td></tr><tr><td>ΔT/Δlat low alt [°C lat<sup>-1</sup>]</td><td>-0.66</td><td>1.14</td></tr><tr><td>ΔT/Δlat high alt [°C lat<sup>-1</sup>]</td><td>-0.92</td><td>0.73</td></tr><tr><td>ΔT/∂z low alt [°C km<sup>-1</sup>]</td><td>5.9</td><td>6.5</td></tr><tr><td>ΔT/Δz high alt [°C km<sup>-1</sup>]</td><td>8.2</td><td>9.8</td></tr><tr><td>MBE [°C] (Obs)</td><td>0.96</td><td>-1.1e-1</td></tr><tr><td>MBE [°C] (MAR)</td><td>4.2e-3</td><td>0.86</td></tr><tr><td>RMSE [°C] (Obs)</td><td>2.8</td><td>4.5</td></tr><tr><td>RMSE [°C] (MAR)</td><td>4.2</td><td>3.8</td></tr></tbody></table>

## Parameterization of mean logarithm accumulation ln(Ā)

To estimate the logarithm of accumulation (ln(Ā)) in Eqs. A1 and A2, we define a linear relationship with the mean annual surface temperature, according to the Clausius Clapeyron law.

In Antarctica (AIS), accumulation records also include surface temperature measurements. These are used to calibrate the dependence of ln(Ā) on surface temperature. For Greenland, no such temperature data are available. Nevertheless, the surface elevation and latitude at the observation sites allow us to compute the parameterized mean surface temperature, as defined earlier.

We derive the following expression to parameterize ln(Ā) in function of the mean annual surface temperature T<sub>s</sub>:

$$ \ln(\bar{A}) = \max \left( \frac{\Delta \ln(\bar{A})}{\Delta T_s} \cdot (T_s - T_{s\_\text{ref}}) + \ln(\bar{A}_{\text{ref}}), \ln(\bar{A}_{\text{min}}) \right) \quad (A8) $$
<figure><img src="image_49.png" /></figure>

The BFGS optimization algorithm is used to determine the optimal rate of variation of the mean logarithm accumulation with respect to the parameterized mean surface temperature. This optimization process aims to minimize the RMSE between the observed and parameterized mean logarithm of accumulation (Fig. A2). The optimized parameters used in the logarithm of the accumulation parameterization can be found in Table A2.

<table><caption>**Table A2.** Parameters used in the logarithm of the accumulation parameterization</caption><thead><tr><th>Parameter</th><th>Greenland</th><th>Antarctica</th></tr></thead><tbody><tr><td>log(A<sub>ref</sub>) [ln(kg<sup>-2</sup>yr<sup>-1</sup>)]</td><td>3.6</td><td>3.4</td></tr><tr><td>T<sub>s_ref</sub> [°C]</td><td>-43</td><td>-59</td></tr><tr><td>ln(A<sub>min</sub>) [ln(kg<sup>-2</sup>yr<sup>-1</sup>)]</td><td>5.09</td><td>3.4</td></tr><tr><td>Δln(A) / ΔT<sub>s</sub> [ln(kg<sup>-2</sup>yr<sup>-1</sup>C<sup>-1</sup>)]</td><td>0.06</td><td>0.06</td></tr><tr><td>MBE [ln(kg<sup>-2</sup>yr<sup>-1</sup>) (Obs)]</td><td>-5.5e-7</td><td>9.3e-3</td></tr><tr><td>MBE [ln(kg<sup>-2</sup>yr<sup>-1</sup>) (MAR)]</td><td>-3.6e-1</td><td>-2.2e-2</td></tr><tr><td>RMSE [ln(kg<sup>-2</sup>yr<sup>-1</sup>) (Obs)]</td><td>0.6</td><td>0.56</td></tr><tr><td>RMSE [ln(kg<sup>-2</sup>yr<sup>-1</sup>) (MAR)]</td><td>0.84</td><td>0.59</td></tr></tbody></table>

## Parameterization of surface density $\rho_s$

In this study, for the initialization, we choose to define an empirical linear relationship between the observed surface density and the mean annual surface temperature, derived from the parameterization based on latitude and elevation.

$$ \rho_s(T_s) = \min \left( \max \left( \left( \frac{\rho_2 - \rho_1}{T_2 - T_1} \right) (T_s - T_1) + \rho_1, \rho_{s\_min} \right), \rho_{s\_max} \right) \quad (A9) $$

The empirical relationship between observed surface density and parameterized surface temperature is then optimized by the BFGS algorithm, separately for Greenland and Antarctica, in order to minimize the RMSE between observed and parameterized surface density. The parameters used in the surface snow density parameterization can be found in Table A3.

<table><caption>**Table A3.** Parameters used in the snow density parameterization</caption><thead><tr><th>Parameter</th><th>Greenland</th><th>Antarctica</th></tr></thead><tbody><tr><td>T<sub>1</sub> [°C]</td><td>-29.0</td><td>-50.2</td></tr><tr><td>T<sub>2</sub> [°C]</td><td>12.0</td><td>-18.0</td></tr><tr><td>ρ<sub>1</sub> [kg m<sup>-3</sup>]</td><td>310.4</td><td>342.0</td></tr><tr><td>ρ<sub>2</sub> [kg m<sup>-3</sup>]</td><td>415.7</td><td>419.3</td></tr><tr><td>ρ<sub>s_min</sub> [kg m<sup>-3</sup>]</td><td>301.5</td><td>341.6</td></tr><tr><td>ρ<sub>s_max</sub> [kg m<sup>-3</sup>]</td><td>398.3</td><td>415.1</td></tr><tr><td>MBE [kg m<sup>-3</sup>] (Obs)</td><td>-9.3e-6</td><td>-6.8e-2</td></tr><tr><td>RMSE [kg m<sup>-3</sup>] (Obs)</td><td>36.36</td><td>29.42</td></tr></tbody></table>

## Parameterization of the snowpack thickness

Our initialization routine allows us to create a snowpack of a uniform thickness with consistent dry snow density properties for the GrIS and AIS. Initializing snowpack thickness can have an influence on the amount of modeled snow melt volumes. A
<figure><img src="image_50.png" /></figure>

uniform snowpack thickness across all regions can introduce biases, especially in areas where the snowpack is naturally thin due to recurrent melt events. This is a concern for the ablation zones of the GrIS. We introduce a variable snowpack thickness in order to better represent these ablation zones of the GrIS. As ablation areas is overall not a concern in Antarctica, this development only applies to the GrIS.

We base our approach on observed density profiles, specifically looking at the depth where density of 600 kg m⁻³ ($dz_{600}$) is reached. This threshold is arbitrary and serves only to differentiate ablation zones from dry snow areas. Ablation zones are characterized by melting and refreezing processes, leading to a thinner snowpack that densifies quickly in comparison to dry snow locations. Thus, with this criteria $dz_{600}$ we can differentiate the ablation zones from the dry snow zones, where densification is slower and less affected by melt-refreeze cycles. After calculating the indicator $dz_{600}$ from density profile trends (to minimize observational noise), we extract the parameterized temperature at the observed points and define the parameterization for snowpack thickness as a function of this variable (Fig. A1).

<figure><img src="image_51.png" /><figcaption>**Figure A1.** Observed snowpack thickness at which 600 kg m⁻³ is reached ($z_{600}$) as a function of the parameterized mean surface temperature and surface elevation, in Greenland. The grey line represents the parameterization of mean snowpack thickness in function the parameterized mean surface temperature.</figcaption></figure>

$$z_{\text{snow}}(T_s) = \min \left( \max \left( z_{\text{snow2}} - \left( \frac{z_{\text{snow2}} - z_{\text{snow1}}}{T_{s2} - T_{s1}} \right) \cdot (T_s - T_{s1}), z_{\text{snow\_min}} \right), z_{\text{snow\_max}} \right) \quad (\text{A10})$$

This results in a parameterization that allows us to distinguish at the initialization the snowpack thickness of the ablation zone from the rest of the GrIS. The parameters used for the snowpack thickness over Greenland are presented in Table A4.
<figure><img src="image_52.png" /><figcaption>**Table A4.** Parameters for snow thickness (Greenland)</figcaption></figure>

<table><thead><tr><td>Parameters</td><td>Greenland</td></tr></thead><tbody><tr><td>T<sub>s1</sub> [°C]</td><td>-14.7</td></tr><tr><td>T<sub>s2</sub> [°C]</td><td>-8.7</td></tr><tr><td>z<sub>snow1</sub> [m]</td><td>0.5</td></tr><tr><td>z<sub>snow2</sub> [m]</td><td>10.0</td></tr><tr><td>z<sub>snow_max</sub> [m]</td><td>0.5</td></tr><tr><td>z<sub>snow_max</sub> [m]</td><td>10.0</td></tr></tbody></table>

<figure><img src="image_53.png" /><figcaption>**Figure A2.** Observed 10 m snow temperature in function of surface elevation and latitude in Greenland (a) and Antarctica (d). Colored lines represent the parameterization of mean annual surface temperature in function of latitude and elevation. Logarithm accumulation observations in function of the parameterized mean annual surface temperature and surface elevation in Greenland (b) and Antarctica (e). Grey line shows the final parameterization in function of the parameterized mean annual surface temperature. Surface density observations in function of the parameterized mean surface temperature and surface elevation in Greenland (c) and Antarctica (f). Grey line shows the parameterized surface density in function of the parameterized mean annual surface temperature.</figcaption></figure>
<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /><figcaption>**Figure A3.** Comparison of observed surface snow density in Antarctica (circle) and Greenland (diamond) as a function of MAR 40 year (1980-2019) average air temperature precipitation weighted (x-axis) and wind speed precipitation weighted (y-axis).</figcaption></figure>
<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /><figcaption>**Figure A4.** Marginal (diagonal) and joint (off-diagonal) distributions of the 30 sets of viscosity parameters ($\eta_0$, $b_\eta$, and $a_\eta$) of the fifth History Matching wave.</figcaption></figure>
<figure><img src="image_58.png" /></figure>

(http://www.cecill.info/index.en.html, last access: 5 September 2023). The ORCHIDEE model code is written in Fortran 90 and is maintained and developed under a Subversion (SVN) version control system at the Institut Pierre-Simon Laplace (IPSL) in France. The ORCHIDAS data assimilation system (in Python) is available at the following address (https://orchidas.lsce.ipsl.fr, last access: 5 November 2025). The python scripts written to generate the analyses and figures for this study are available at https://doi.org/10.5281/zenodo.17569822. The datasets and model outputs used for this study are available at https://doi.org/10.5281/zenodo.17608287. The outputs from the HM calibration procedure are available at https://doi.org/10.5281/zenodo.17523380. The MAR outputs used to force the snowpack model are available at http://ftp.climato.be/fettweis/MARv3.13/Greenland/ (last access: 5 November 2025) and at http://ftp.climato.be/fettweis/MARv3.12/Antarctica/ (last access: 5 November 2025) for the Greenland and Antarctic ice sheets respectively. The SUMUP dataset used for this study is available online at https://doi.org/10.18739/A2M61BR5M. All the python and fortran scripts for the initialization procedure are publicly available via: https://gitlab.in2p3.fr/ipsl/projets/awaca/models/init-snowpack.git

**Author contributions.** PC, SC and CA co-designed the overall work. CA developed the early phases of the initialization procedure for the AIS, which was further adapted to the Greenland ice sheet by PC. PC implemented the presented developments in ORCHIDEE. PC conducted the simulations, processed the data, analyzed the model outputs, with the help of SC and CA. CD implemented the new refreezing scheme in the snowpack model. SB and NR provided expertise on the History Matching (HM) calibration procedure. PC performed the calibration of snow density with HM, with the advice of SB. PC wrote the manuscript, with contributions from SC and CA, and produced the figures. All authors provided comments on the manuscript.

**Competing interests.** The contact author has declared that none of the authors had any competing interests.

**Acknowledgements.** This work benefited from the computing facilities and post-processing softwares from the IPSL Data and Computing Center ESPRI, which is supported by CNRS, Sorbonne University, CNES and Ecole Polytechnique. We would like to thank Xavier Fettweis for providing the MAR outputs for the Greenland and Antarctic ice sheets and Laurent Arnaud for providing the density profiles from Dome C. We also thank Marie Dumont and Vincent Favier, as well as Charles Amory and all the participants involved in the SNOW working group, composed of members of the Institut Pierre Simon Laplace (IPSL, France) and the Institut des Géosciences de l'Environnement (IGE, France), for their numerous and fruitful discussions. We would like to thank the core ORCHIDEE team for maintaining the model, and especially Josefine Ghattas for helping merge the new developments of the snowpack scheme into the trunk version of the ORCHIDEE model. We are grateful to the ORCHIDAS team for maintaining the software in which History Matching is implemented. Nina Raoult is funded by the Destination Earth initiative of the European Commission. Simon Beylat has been supported by a scholarship from CNRS under the Melbourne-CNRS joint doctoral programme.

**Financial support.** This study has received funding from the Agence Nationale de la Recherche - France 2030 as part of the PEPR TRACCS program (grant nos. ANR-22-EXTR-0010 and ANR-22-EXTR-0008).
<figure><img src="image_59.png" /></figure>

References

Agosta, C., Amory, C., Kittel, C., Orsi, A., Favier, V., Gallée, H., van den Broeke, M. R., Lenaerts, J. T. M., van Wessem, J. M., van de Berg, W. J., and Fettweis, X.: Estimation of the Antarctic surface mass balance using the regional climate model MAR (1979–2015) and identification of dominant processes, The Cryosphere, 13, 281–296, https://doi.org/10.5194/tc-13-281-2019, publisher: Copernicus GmbH, 2019.

Alley, R.: Firn densification by grain-boundary sliding: a first model, https://hal.science/file/index/docid/226281/filename/ajp-jphyscol198748C135.pdf, 1987.

Alley, R. B., Clark, P. U., Huybrechts, P., and Joughin, I.: Ice-Sheet and Sea-Level Changes, Science, 310, 456–460, https://doi.org/10.1126/science.1114613, publisher: American Association for the Advancement of Science, 2005.

Amory, C., Kittel, C., Le Toumelin, L., Agosta, C., Delhasse, A., Favier, V., and Fettweis, X.: Performance of MAR (v3.11) in simulating the drifting-snow climate and surface mass balance of Adélie Land, East Antarctica, Geoscientific Model Development, 14, 3487–3510, https://doi.org/10.5194/gmd-14-3487-2021, publisher: Copernicus GmbH, 2021.

Amory, C., Buizert, C., Buzzard, S., Case, E., Clerx, N., Culberg, R., Datta, R. T., Dey, R., Drews, R., Dunmire, D., Eayrs, C., Hansen, N., Humbert, A., Kaitheri, A., Keegan, K., Kuipers Munneke, P., Lenaerts, J. T. M., Lhermitte, S., Mair, D., McDowell, I., Mejia, J., Meyer, C. R., Morris, E., Moser, D., Oraschewski, F. M., Pearce, E., de Roda Husman, S., Schlegel, N.-J., Schultz, T., Simonsen, S. B., Stevens, C. M., Thomas, E. R., Thompson-Munson, M., Wever, N., Wouters, B., and The Firn Symposium team: Firn on ice sheets, Nature Reviews Earth & Environment, 5, 79–99, https://doi.org/10.1038/s43017-023-00507-9, publisher: Nature Publishing Group, 2024.

Anderson, E. A.: A point energy and mass balance model of a snow cover, Ph.D. thesis, Stanford University, USA, https://repository.library.noaa.gov/view/noaa/6392/, 1976.

Arnaud, L., Lipenkov, V., Barnola, J. M., Gay, M., and Duval, P.: Modelling of the densification of polar firn: characterization of the snow-firn transition, Annals of Glaciology, 26, 39–44, https://doi.org/10.3189/1998AoG26-1-39-44, 1998.

Arthern, R. J., Vaughan, D. G., Rankin, A. M., Mulvaney, R., and Thomas, E. R.: In situ measurements of Antarctic snow compaction compared with predictions of models, Journal of Geophysical Research: Earth Surface, 115, https://doi.org/10.1029/2009JF001306, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2009JF001306, 2010.

Baker, E., Harper, A. B., Williamson, D., and Challenor, P.: Emulation of high-resolution land surface models using sparse Gaussian processes with application to JULES, Geoscientific Model Development, 15, 1913–1929, https://doi.org/10.5194/gmd-15-1913-2022, publisher: Copernicus GmbH, 2022.

Barnola, J. M., Pimienta, P., Raynaud, D., and Korotkevich, Y. S.: CO2-climate relationship as deduced from the Vostok ice core: a re-examination based on new measurements and on a re-evaluation of the air dating, Tellus B, 43, 83–90, https://doi.org/10.1034/j.1600-0889.1991.t01-1-00002.x, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1034/j.1600-0889.1991.t01-1-00002.x, 1991.

Barrere, M., Domine, F., Decharme, B., Morin, S., Vionnet, V., and Lafaysse, M.: Evaluating the performance of coupled snow–soil models in SURFEXv8 to simulate the permafrost thermal regime at a high Arctic site, Geoscientific Model Development, 10, 3461, https://doi.org/10.5194/gmd-10-3461-2017, 2017.

Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning: Part I: numerical model, Cold Regions Science and Technology, 35, 123–145, https://doi.org/10.1016/S0165-232X(02)00074-5, 2002.
<figure><img src="image_60.png" /></figure>

Boone, A. and Etchevers, P.: An Intercomparison of Three Snow Schemes of Varying Complexity Coupled to the Same Land Surface Model: Local-Scale Evaluation at an Alpine Site, Journal of Hydrometeorology, 2, 374–394, https://doi.org/10.1175/1525-7541(2001)002<0374:AIOTSS>2.0.CO;2, publisher: American Meteorological Society Section: Journal of Hydrometeorology, 2001.

Boucher, O., Servonnat, J., Albright, A. L., Aumont, O., Balkanski, Y., Bastrikov, V., Bekki, S., Bonnet, R., Bony, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Caubel, A., Cheruy, F., Codron, F., Cozic, A., Cugnet, D., D'Andrea, F., Davini, P., de Lavergne, C., Denvil, S., Deshayes, J., Devilliers, M., Ducharne, A., Dufresne, J.-L., Dupont, E., Éthé, C., Fairhead, L., Falletti, L., Flavoni, S., Foujols, M.-A., Gardoll, S., Gastineau, G., Ghattas, J., Grandpeix, J.-Y., Guenet, B., Guez, Lionel, E., Guilyardi, E., Guimberteau, M., Hauglustaine, D., Hourdin, F., Idelkadi, A., Joussaume, S., Kageyama, M., Khodri, M., Krinner, G., Lebas, N., Levavasseur, G., Lévy, C., Li, L., Lott, F., Lurton, T., Luyssaert, S., Madec, G., Madeleine, J.-B., Maignan, F., Marchand, M., Marti, O., Mellul, L., Meurdesoif, Y., Mignot, J., Musat, I., Ottlé, C., Peylin, P., Planton, Y., Polcher, J., Rio, C., Rochetin, N., Rousset, C., Sepulchre, P., Sima, A., Swingedouw, D., Thiéblemont, R., Traore, A. K., Vancoppenolle, M., Vial, J., Vialard, J., Viovy, N., and Vuichard, N.: Presentation and Evaluation of the IPSL-CM6A-LR Climate Model, Journal of Advances in Modeling Earth Systems, 12, e2019MS002010, https://doi.org/10.1029/2019MS002010, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2019MS002010, 2020.

Box, J. E., Hubbard, A., Bahr, D. B., Colgan, W. T., Fettweis, X., Mankoff, K. D., Wehrlé, A., Noël, B., van den Broeke, M. R., Wouters, B., Bjørk, A. A., and Fausto, R. S.: Greenland ice sheet climate disequilibrium and committed sea-level rise, Nature Climate Change, 12, 808–813, https://doi.org/10.1038/s41558-022-01441-2, publisher: Nature Publishing Group, 2022.

Brun, E., Martin, Simon, V., Gendre, C., and Coleou, C.: An Energy and Mass Model of Snow Cover Suitable for Operational Avalanche Forecasting, Journal of Glaciology, 35, 333–342, https://doi.org/10.3189/S0022143000009254, 1989.

Buizert, C., Fudge, T. J., Roberts, W. H. G., Steig, E. J., Sherriff-Tadano, S., Ritz, C., Lefebvre, E., Edwards, J., Kawamura, K., Oyabu, I., Motoyama, H., Kahle, E. C., Jones, T. R., Abe-Ouchi, A., Obase, T., Martin, C., Corr, H., Severinghaus, J. P., Beaudette, R., Epifanio, J. A., Brook, E. J., Martin, K., Chappellaz, J., Aoki, S., Nakazawa, T., Sowers, T. A., Alley, R. B., Ahn, J., Sigl, M., Severi, M., Dunbar, N. W., Svensson, A., Fegyveresi, J. M., He, C., Liu, Z., Zhu, J., Otto-Bliesner, B. L., Lipenkov, V. Y., Kageyama, M., and Schwander, J.: Antarctic surface temperature and elevation during the Last Glacial Maximum, Science, 372, 1097–1101, https://doi.org/10.1126/science.abd2897, publisher: American Association for the Advancement of Science, 2021.

Chalita, S. and Le Treut, H.: The albedo of temperate and boreal forest and the Northern Hemisphere climate: a sensitivity experiment using the LMD GCM, Climate Dynamics, 10, 231–240, https://doi.org/10.1007/BF00208990, 1994.

Charbit, S., Dumas, C., Maignan, F., Ottlé, C., Raoult, N., Fettweis, X., and Conesa, P.: Modelling snowpack on ice surfaces with the ORCHIDEE land surface model: application to the Greenland ice sheet, The Cryosphere, 18, 5067–5099, https://doi.org/10.5194/tc-18-5067-2024, publisher: Copernicus GmbH, 2024.

Cheruy, F., Ducharne, A., Hourdin, F., Musat, I., Vignon, Gastineau, G., Bastrikov, V., Vuichard, N., Diallo, B., Dufresne, J.-L., Ghattas, J., Grandpeix, J.-Y., Idelkadi, A., Mellul, L., Maignan, F., Ménégoz, M., Ottlé, C., Peylin, P., Servonnat, J., Wang, F., and Zhao, Y.: Improved Near-Surface Continental Climate in IPSL-CM6A-LR by Combined Evolutions of Atmospheric and Land Surface Physics, Journal of Advances in Modeling Earth Systems, 12, e2019MS002005, https://doi.org/10.1029/2019MS002005, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2019MS002005, 2020.

Coléou, C. and Lesaffre, B.: Irreducible water saturation in snow: experimental results in a cold laboratory, Annals of Glaciology, 26, 64–68, https://doi.org/10.3189/1998AoG26-1-64-68, 1998.

Couvreux, F., Hourdin, F., Williamson, D., Roehrig, R., Volodina, V., Villefranque, N., Rio, C., Audouin, O., Salter, J., Bazile, E., Brient, F., Favot, F., Honnert, R., Lefebvre, M.-P., Madeleine, J.-B., Rodier, Q., and Xu, W.: Process-Based Climate Model Development Har-
<figure><img src="image_61.png" /></figure>

765

nessing Machine Learning: I. A Calibration Tool for Parameterization Improvement, Journal of Advances in Modeling Earth Systems, 13, e2020MS002217, https://doi.org/10.1029/2020MS002217, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2020MS002217, 2021.

de la Peña, S., Howat, I. M., Nienow, P. W., van den Broeke, M. R., Mosley-Thompson, E., Price, S. F., Mair, D., Noël, B., and Sole, A. J.: Changes in the firn structure of the western Greenland Ice Sheet caused by recent warming, The Cryosphere, 9, 1203–1211, https://doi.org/10.5194/tc-9-1203-2015, publisher: Copernicus GmbH, 2015.

770

Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P., and Morin, S.: Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model, The Cryosphere, 10, 853–877, https://doi.org/10.5194/tc-10-853-2016, publisher: Copernicus GmbH, 2016.

Domine, F., Barrere, M., and Morin, S.: The growth of shrubs on high Arctic tundra at Bylot Island: impact on snow physical properties and permafrost thermal regime, Biogeosciences, 13, 6471–6486, https://doi.org/10.5194/bg-13-6471-2016, publisher: Copernicus GmbH, 2016.

775

Ekaykin, A. A., Bolshunov, A. V., Lipenkov, V. Y., Scheinert, M., Eberlein, L., Brovkov, E., Popov, S. V., and Turkeev, A. V.: First glacio-logical investigations at Ridge B, central East Antarctica, Antarctic Science, 33, 418–427, https://doi.org/10.1017/S0954102021000171, 2021.

Ekaykin, A. A., Tchikhatchev, K. B., Veres, A. N., Lipenkov, V. Y., Tebenkova, N. A., and Turkeev, A. V.: Vertical profile of snow-firn density in the vicinity of Vostok station, Central Antarctica, ë, 62, 504–511, https://doi.org/10.31857/s2076673422040147, publisher: Nauka, 2023.

780

Fausto, R. S., Box, J. E., Vandegruch, B., van As, D., Steffen, K., MacFerrin, M. J., Machguth, H., Colgan, W., Koenig, L. S., McGrath, D., Charalampidis, C., and Braithwaite, R. J.: A Snow Density Dataset for Improving Surface Boundary Conditions in Greenland Ice Sheet Firn Modeling, Frontiers in Earth Science, 6, https://doi.org/10.3389/feart.2018.00051, publisher: Frontiers, 2018.

Favier, V., Agosta, C., Parouty, S., Durand, G., Delaygue, G., Gallée, H., Drouet, A.-S., Trouvilliez, A., and Krinner, G.: An updated and quality controlled surface mass balance dataset for Antarctica, The Cryosphere, 7, 583–597, https://doi.org/10.5194/tc-7-583-2013, publisher: Copernicus GmbH, 2013.

785

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033, https://doi.org/10.5194/tc-11-1015-2017, publisher: Copernicus GmbH, 2017.

790

Fettweis, X., Ambroise, B., David, P.-M., Ghilain, N., Paul, P., and Wuest, C.: Évolution actuelle (1960-2021) de l'enneigement dans les Vosges à l'aide du modèle régional du climat MAR, Bulletin de la Société Géographique de Liège, https://doi.org/10.25518/0770-7576.7049, 2023.

Feulner, G., Rahmstorf, S., Levermann, A., and Volkwardt, S.: On the Origin of the Surface Air Temperature Difference between the Hemispheres in Earth's Present-Day Climate, Journal of Climate, 26, 7136–7150, https://doi.org/10.1175/JCLI-D-12-00636.1, publisher: American Meteorological Society Section: Journal of Climate, 2013.

795

Forsberg, R., Sørensen, L., and Simonsen, S.: Greenland and Antarctica Ice Sheet Mass Changes and Effects on Global Sea Level, Surveys in Geophysics, 38, 89–104, https://doi.org/10.1007/s10712-016-9398-7, 2017.

Genthon, C., Six, D., Scarchilli, C., Ciardini, V., and Frezzotti, M.: Meteorological and snow accumulation gradients across Dome C, East Antarctic plateau, International Journal of Climatology, 36, 455–466, https://doi.org/10.1002/joc.4362, _eprint: https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/joc.4362, 2016.

800
<figure><img src="image_62.png" /></figure>

Goujon, C., Barnola, J.-M., and Ritz, C.: Modeling the densification of polar firn including heat diffusion: Application to close-off char- acteristics and gas isotopic fractionation for Antarctica and Greenland sites, Journal of Geophysical Research: Atmospheres, 108, https://doi.org/10.1029/2002JD003319, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2002JD003319, 2003.

805
Hanna, E., Huybrechts, P., Janssens, I., Cappelen, J., Steffen, K., and Stephens, A.: Runoff and mass balance of the Green-
land ice sheet: 1958–2003, Journal of Geophysical Research: Atmospheres, 110, https://doi.org/10.1029/2004JD005641, _eprint:
https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2004JD005641, 2005.

Helsen, M. M., van de Wal, R. S. W., Reerink, T. J., Bintanja, R., Madsen, M. S., Yang, S., Li, Q., and Zhang, Q.: On the impor-
tance of the albedo parameterization for the mass balance of the Greenland ice sheet in EC-Earth, The Cryosphere, 11, 1949–1965,
https://doi.org/10.5194/tc-11-1949-2017, publisher: Copernicus GmbH, 2017.

Herron, M. M. and Jr, C. C. L.: Firn Densification: An Empirical Model, Journal of Glaciology, 25, 373–385,
https://doi.org/10.3189/S0022143000015239, 1980.

1815
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Sim-
mons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren,
P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J.,
Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Vil-
laume, S., and Thépaut, J.-N.: The ERA5 global reanalysis, Quarterly Journal of the Royal Meteorological Society, 146, 1999–2049,
https://doi.org/10.1002/qj.3803, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/qj.3803, 2020.

820
Hourdin, F., Williamson, D., Rio, C., Couvreux, F., Roehrig, R., Villefranque, N., Musat, I., Fairhead, L., Diallo, F. B., and Volo-
dina, V.: Process-Based Climate Model Development Harnessing Machine Learning: II. Model Calibration From Single Column
to Global, Journal of Advances in Modeling Earth Systems, 13, e2020MS002225, https://doi.org/10.1029/2020MS002225, _eprint:
https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020MS002225, 2021.

820

Hourdin, F., Ferster, B., Deshayes, J., Mignot, J., Musat, I., and Williamson, D.: Toward machine-assisted tuning avoiding the underestimation of uncertainty in climate change projections, Science Advances, 9, eadf2758, https://doi.org/10.1126/sciadv.adf2758, publisher: American Association for the Advancement of Science, 2023.

825

Howat, I. M.: Temporal variability in snow accumulation and density at Summit Camp, Greenland ice sheet, Journal of Glaciology, 68, 1076–1084, https://doi.org/10.1017/jog.2022.21, 2022.

Keenan, E., Wever, N., Dattler, M., Lenaerts, J. T. M., Medley, B., Kuipers Munneke, P., and Reijmer, C.: Physics-based SNOWPACK model improves representation of near-surface Antarctic snow and firn density, The Cryosphere, 15, 1065–1085, https://doi.org/10.5194/tc-15-1065-2021, publisher: Copernicus GmbH, 2021.

830

Kipfstuhl, S., Faria, S. H., Azuma, N., Freitag, J., Hamann, I., Kaufmann, P., Miller, H., Weiler, K., and Wilhelms, F.: Evidence of dy-
namic recrystallization in polar firn, Journal of Geophysical Research: Solid Earth, 114, https://doi.org/10.1029/2008JB005583, _eprint:
https://onlinelibrary.wiley.com/doi/pdf/10.1029/2008JB005583, 2009.

Kittel, C., Amory, C., Agosta, C., Jourdain, N. C., Hofer, S., Delhasse, A., Doutreloup, S., Huot, P.-V., Lang, C., Fichefet, T., and Fettweis, X.: Diverging future surface mass balance between the Antarctic ice shelves and grounded ice sheet, The Cryosphere, 15, 1215–1236, https://doi.org/10.5194/tc-15-1215-2021, publisher: Copernicus GmbH, 2021.

835

Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S., and Prentice, I. C.: A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system, Global Biogeochemical Cycles, 19, https://doi.org/10.1029/2003GB002199, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2003GB002199, 2005.
<figure><img src="image_63.png" /></figure>

840 Kuipers Munneke, P., Reijmer, C. H., van den Broeke, M. R., König-Langlo, G., Stammes, P., and Knap, W. H.: Analysis of clear-sky Antarctic snow albedo using observations and radiative transfer modeling, Journal of Geophysical Research: Atmospheres, 113, https://doi.org/10.1029/2007JD009653, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2007JD009653, 2008.

Lafaysse, M., Cluzet, B., Dumont, M., Lejeune, Y., Vionnet, V., and Morin, S.: A multiphysical ensemble system of numerical snow modeling, The Cryosphere, 11, 1173–1198, https://doi.org/10.5194/tc-11-1173-2017, publisher: Copernicus GmbH, 2017.

845 Lambin, C., Fettweis, X., Kittel, C., Fonder, M., and Ernst, D.: Assessment of future wind speed and wind power changes over South Greenland using the MAR regional climate model, International Journal of Climatology, https://doi.org/10.1002/joc.7795, publisher: John Wiley & Sons, Hoboken, United States - New Jersey, 2022.

Leduc-Leballeur, M., Picard, G., Mialon, A., Arnaud, L., Lefebvre, E., Possenti, P., and Kerr, Y.: Modeling L-Band Brightness Temperature at Dome C in Antarctica and Comparison With SMOS Observations, IEEE Transactions on Geoscience and Remote Sensing, 53, 4022–4032, https://doi.org/10.1109/TGRS.2015.2388790, 2015.

Leduc-Leballeur, M., Picard, G., Macelloni, G., Mialon, A., and Kerr, Y. H.: Melt in Antarctica derived from Soil Moisture and Ocean Salinity (SMOS) observations at L band, The Cryosphere, 14, 539–548, https://doi.org/10.5194/tc-14-539-2020, publisher: Copernicus GmbH, 2020.

Lenaerts, J. T. M., van den Broeke, M. R., Déry, S. J., van Meijgaard, E., van de Berg, W. J., Palm, S. P., and Sanz Rodrigo, J.: Modeling drifting snow in Antarctica with a regional climate model: 1. Methods and model evaluation, Journal of Geophysical Research: Atmospheres, 117, https://doi.org/10.1029/2011JD016145, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2011JD016145, 2012.

Ligtenberg, S. R. M., Helsen, M. M., and van den Broeke, M. R.: An improved semi-empirical model for the densification of Antarctic firn, The Cryosphere, 5, 809–819, https://doi.org/10.5194/tc-5-809-2011, publisher: Copernicus GmbH, 2011.

Liston, G. E. and Sturm, M.: Winter Precipitation Patterns in Arctic Alaska Determined from a Blowing-Snow Model and Snow-Depth Observations, Journal of Hydrometeorology, 3, 646–659, https://doi.org/10.1175/1525-7541(2002)003<0646:WPPIAA>2.0.CO;2, publisher: American Meteorological Society Section: Journal of Hydrometeorology, 2002.

Loeppky, J. L., Sacks, J., and Welch, W. J.: Choosing the Sample Size of a Computer Experiment: A Practical Guide, Technometrics, 51, 366–376, https://doi.org/10.1198/TECH.2009.08040, publisher: ASA Website _eprint: https://doi.org/10.1198/TECH.2009.08040, 2009.

MacBean, N., Bacour, C., Raoult, N., Bastrikov, V., Koffi, E. N., Kuppel, S., Maignan, F., Ottlé, C., Peaucelle, M., Santaren, D., and Peylin, P.: Quantifying and Reducing Uncertainty in Global Carbon Cycle Predictions: Lessons and Perspectives From 15 Years of Data Assimilation Studies With the ORCHIDEE Terrestrial Biosphere Model, Global Biogeochemical Cycles, 36, e2021GB007177, https://doi.org/10.1029/2021GB007177, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2021GB007177, 2022.

Marquetto, L., Kaspari, S., and Cardia Simões, J.: Refractory black carbon (rBC) variability in a 47-year West Antarctic snow and firn core, The Cryosphere, 14, 1537–1554, https://doi.org/10.5194/tc-14-1537-2020, publisher: Copernicus GmbH, 2020.

McNeall, D., Robertson, E., and Wiltshire, A.: Constraining the carbon cycle in JULES-ES-1.0, Geoscientific Model Development, 17, 1059–1089, https://doi.org/10.5194/gmd-17-1059-2024, publisher: Copernicus GmbH, 2024.

McNeall, D. J., Challenor, P. G., Gattiker, J. R., and Stone, E. J.: The potential of an observational data set for calibration of a computationally expensive computer model, Geoscientific Model Development, 6, 1715–1728, https://doi.org/10.5194/gmd-6-1715-2013, publisher: Copernicus GmbH, 2013.

Morris, M. D. and Mitchell, T. J.: Exploratory designs for computational experiments, Journal of Statistical Planning and Inference, 43, 381–402, https://doi.org/10.1016/0378-3758(94)00035-T, 1995.
<figure><img src="image_64.png" /></figure>

Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.: Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, Proceedings of the National Academy of Sciences, 116, 9239–9244, https://doi.org/10.1073/pnas.1904242116, publisher: Proceedings of the National Academy of Sciences, 2019.

880 Mousavi, M., Colliander, A., Miller, J. Z., Entekhabi, D., Johnson, J. T., Shuman, C. A., Kimball, J. S., and Courville, Z. R.: Evaluation of Surface Melt on the Greenland Ice Sheet Using SMAP L-Band Microwave Radiometry, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 14, 11439–11449, https://doi.org/10.1109/JSTARS.2021.3124229, 2021.

885 Oppenheimer, M., Glavovic, and Hinkel: Sea Level Rise and Implications for Low-Lying Islands, Coasts and Communities, The Ocean and Cryosphere in a Changing Climate: Special Report of the Intergovernmental Panel on Climate Change, pp. 321–446, https://doi.org/10.1017/9781009157964.006, 2022.

890 Otosaka, I. N., Shepherd, A., Ivins, E. R., Schlegel, N.-J., Amory, C., van den Broeke, M. R., Horwath, M., Joughin, I., King, M. D., Krinner, G., Nowicki, S., Payne, A. J., Rignot, E., Scambos, T., Simon, K. M., Smith, B. E., Sørensen, L. S., Velicogna, I., Whitehouse, P. L., A., G., Agosta, C., Ahlstrøm, A. P., Blazquez, A., Colgan, W., Engdahl, M. E., Fettweis, X., Forsberg, R., Gallée, H., Gardner, A., Gilbert, L., Gourmelen, N., Groh, A., Gunter, B. C., Harig, C., Helm, V., Khan, S. A., Kittel, C., Konrad, H., Langen, P. L., Lecavalier, B. S., Liang, C.-C., Loomis, B. D., McMillan, M., Melini, D., Mernild, S. H., Mottram, R., Mouginot, J., Nilsson, J., Noël, B., Pattle, M. E., Peltier, W. R., Pie, N., Roca, M., Sasgen, I., Save, H. V., Seo, K.-W., Scheuchl, B., Schrama, E. J. O., Schröder, L., Simonsen, S. B., Slater, T., Spada, G., Sutterley, T. C., Vishwakarma, B. D., van Wessem, J. M., Wiese, D., van der Wal, W., and Wouters, B.: Mass balance of the Greenland and Antarctic ice sheets from 1992 to 2020, Earth System Science Data, 15, 1597–1616, https://doi.org/10.5194/essd-15-1597-2023, publisher: Copernicus GmbH, 2023.

895 Pahaut, E.: La métamorphose des cristaux de neige (Snow crystal metamorphosis), Monographies de la Météorologie Nationale, 1976.

Polashenski, C., Courville, Z., Benson, C., Wagner, A., Chen, J., Wong, G., Hawley, R., and Hall, D.: Observations of pronounced Greenland ice sheet firn warming and implications for runoff production, Geophysical Research Letters, 41, 4238–4246, https://doi.org/10.1002/2014GL059806, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2014GL059806, 2014.

900 Raoult, N., Charbit, S., Dumas, C., Maignan, F., Ottlé, C., and Bastrikov, V.: Improving modelled albedo over the Greenland ice sheet through parameter optimisation and MODIS snow albedo retrievals, The Cryosphere, 17, 2705–2724, https://doi.org/10.5194/tc-17-2705-2023, publisher: Copernicus GmbH, 2023.

Raoult, N., Beylat, S., Salter, J. M., Hourdin, F., Bastrikov, V., Ottlé, C., and Peylin, P.: Exploring the potential of history matching for land surface model calibration, Geoscientific Model Development, 17, 5779–5801, https://doi.org/10.5194/gmd-17-5779-2024, publisher: Copernicus GmbH, 2024.

905 Ridder, K. D. and Schayes, G.: The IAGL Land Surface Model, Journal of Applied Meteorology and Climatology, 36, 167–182, https://doi.org/10.1175/1520-0450(1997)036<0167:TILSM>2.0.CO;2, publisher: American Meteorological Society Section: Journal of Applied Meteorology and Climatology, 1997.

Rignot, E., Mouginot, J., Scheuchl, B., van den Broeke, M., van Wessem, M. J., and Morlighem, M.: Four decades of Antarctic Ice Sheet mass balance from 1979–2017, Proceedings of the National Academy of Sciences, 116, 1095–1103, https://doi.org/10.1073/pnas.1812883116, publisher: Proceedings of the National Academy of Sciences, 2019.

910 Shepherd, A., Ivins, E., Rignot, E., Smith, B., van den Broeke, M., Velicogna, I., Whitehouse, P., Briggs, K., Joughin, I., Krinner, G., Nowicki, S., Payne, T., Scambos, T., Schlegel, N., A, G., Agosta, C., Ahlstrøm, A., Babonis, G., Barletta, V., Blazquez, A., Bonin, J., Csatho, B., Cullather, R., Felikson, D., Fettweis, X., Forsberg, R., Gallee, H., Gardner, A., Gilbert, L., Groh, A., Gunter, B., Hanna, E., Harig, C., Helm, V., Horvath, A., Horwath, M., Khan, S., Kjeldsen, K. K., Konrad, H., Langen, P., Lecavalier, B., Loomis, B., Luthcke, S.,
<figure><img src="image_65.png" /></figure>

915 McMillan, M., Melini, D., Mernild, S., Mohajerani, Y., Moore, P., Mouginot, J., Moyano, G., Muir, A., Nagler, T., Nield, G., Nilsson, J., Noel, B., Otosaka, I., Pattle, M. E., Peltier, W. R., Pie, N., Rietbroek, R., Rott, H., Sandberg-Sørensen, L., Sasgen, I., Save, H., Scheuchl, B., Schrama, E., Schröder, L., Seo, K.-W., Simonsen, S., Slater, T., Spada, G., Sutterley, T., Talpe, M., Tarasov, L., van de Berg, W. J., van der Wal, W., van Wessem, M., Vishwakarma, B. D., Wiese, D., Wouters, B., and The IMBIE team: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 558, 219–222, https://doi.org/10.1038/s41586-018-0179-y, publisher: Nature Publishing Group, 2018.

920 Skiles, S. M., Flanner, M., Cook, J. M., Dumont, M., and Painter, T. H.: Radiative forcing by light-absorbing particles in snow, Nature Climate Change, 8, 964–971, https://doi.org/10.1038/s41558-018-0296-5, publisher: Nature Publishing Group, 2018.

Steger, C. R., Reijmer, C. H., van den Broeke, M. R., Wever, N., Forster, R. R., Koenig, L. S., Kuipers Munneke, P., Lehning, M., Lhermitte, S., Ligtenberg, S. R. M., Miège, C., and Noël, B. P. Y.: Firn Meltwater Retention on the Greenland Ice Sheet: A Model Comparison, Frontiers in Earth Science, 5, https://doi.org/10.3389/feart.2017.00003, publisher: Frontiers, 2017.

925 Stevens, C. M., Verjans, V., Lundin, J. M. D., Kahle, E. C., Horlings, A. N., Horlings, B. I., and Waddington, E. D.: The Community Firn Model (CFM) v1.0, Geoscientific Model Development, 13, 4355–4377, https://doi.org/10.5194/gmd-13-4355-2020, publisher: Copernicus GmbH, 2020.

930 Sturm, M., Glen E. L., Carl S. B., and Holmgren, J.: Characteristics and Growth of a Snowdrift in Arctic Alaska, U.S.A., Arctic, Antarctic, and Alpine Research, 33, 319–329, https://doi.org/10.1080/15230430.2001.12003436, publisher: Taylor & Francis _eprint: https://doi.org/10.1080/15230430.2001.12003436, 2001.

Valt, M., Guyennon, N., Salerno, F., Petrangeli, A. B., Salvatori, R., Cianfarra, P., and Romano, E.: Predicting new snow density in the Italian Alps: A variability analysis based on 10 years of measurements, Hydrological Processes, 32, 3174–3187, https://doi.org/10.1002/hyp.13249, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/hyp.13249, 2018.

935 van Angelen, J. H., M. Lenaerts, J. T., van den Broeke, M. R., Fettweis, X., and van Meijgaard, E.: Rapid loss of firn pore space accelerates 21st century Greenland mass loss, Geophysical Research Letters, 40, 2109–2113, https://doi.org/10.1002/grl.50490, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/grl.50490, 2013.

van den Broeke, M.: Depth and Density of the Antarctic Firn Layer, Arctic, Antarctic, and Alpine Research, 40, 432–438, https://doi.org/10.1657/1523-0430(07-021)[BROEKE]2.0.CO;2, publisher: Taylor & Francis _eprint: https://www.tandfonline.com/doi/pdf/10.1657/1523-0430%2807-021%29%5BBROEKE%5D2.0.CO%3B2, 2008.

940 van den Broeke, M. R., Enderlin, E. M., Howat, I. M., Kuipers Munneke, P., Noël, B. P. Y., van de Berg, W. J., van Meijgaard, E., and Wouters, B.: On the recent contribution of the Greenland ice sheet to sea level change, The Cryosphere, 10, 1933–1946, https://doi.org/10.5194/tc-10-1933-2016, publisher: Copernicus GmbH, 2016.

van Kampenhout, L., Lenaerts, J. T. M., Lipscomb, W. H., Sacks, W. J., Lawrence, D. M., Slater, A. G., and van den Broeke, M. R.: Improving the Representation of Polar Snow and Firn in the Community Earth System Model, Journal of Advances in Modeling Earth Systems, 9, 2583–2600, https://doi.org/10.1002/2017MS000988, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2017MS000988, 2017.

945 Vandegrusk, B., Fausto, R. S., Langen, P. L., van As, D., MacFerrin, M., Colgan, W. T., Ingeman-Nielsen, T., Steffen, K., Jensen, N. S., Møller, M. T., and Box, J. E.: Drivers of Firn Density on the Greenland Ice Sheet Revealed by Weather Station Observations and Modeling, Journal of Geophysical Research: Earth Surface, 123, 2563–2576, https://doi.org/10.1029/2017JF004597, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2017JF004597, 2018.

Vandegrusk, B., MacFerrin, M., Machguth, H., Colgan, W. T., van As, D., Heilig, A., Stevens, C. M., Charalampidis, C., Fausto, R. S., Morris, E. M., Mosley-Thompson, E., Koenig, L., Montgomery, L. N., Miège, C., Simonsen, S. B., Ingeman-Nielsen, T., and Box, J. E.: Firn data
<figure><img src="image_66.png" /></figure>

compilation reveals widespread decrease of firn air content in western Greenland, The Cryosphere, 13, 845–859, https://doi.org/10.5194/tc-13-845-2019, publisher: Copernicus GmbH, 2019.

955 Vandecrux, B., Mottram, R., Langen, P. L., Fausto, R. S., Olesen, M., Stevens, C. M., Verjans, V., Leeson, A., Ligtenberg, S., Kuipers Munneke, P., Marchenko, S., van Pelt, W., Meyer, C. R., Simonsen, S. B., Heilig, A., Samimi, S., Marshall, S., Machguth, H., MacFerrin, M., Niwano, M., Miller, O., Voss, C. I., and Box, J. E.: The firn meltwater Retention Model Intercomparison Project (RetMIP): evaluation of nine firn models at four weather station sites on the Greenland ice sheet, The Cryosphere, 14, 3785–3810, https://doi.org/10.5194/tc-14-3785-2020, publisher: Copernicus GmbH, 2020.

960 Vandecrux, B., Amory, C., Ahlstrøm, A. P., Akers, P. D., Albert, M., Alley, R. B., Arnaud, L., Bales, R., Benson, C., Box, J. E., Buizert, C., Charalampidis, C., Clerx, N., Dibb, J. E., Covi, F., Denis, G., Ding, M., Eisen, O., Fausto, R., Fernandoy, F., Freitag, J., Gerland, S., Harper, J., Hawley, R. L., Hock, R., How, P., Hubbard, B., Humphrey, N., Iizuka, Y., Isaksson, E., Kameda, T., Karlsson, N. B., Kawakami, K., Kjær, H. A., Kuipers_Munneke, P., Lewis, G., MacFerrin, M., Machguth, H., Mankoff, K. D., McConnell, J. R., Medley, B., Morris, E., Mosley-Thompson, E., Mulvaney, R., Niwano, M., Osterberg, E., Otosaka, I., Picard, G., Polashenski, C., Renormalm, A., Rutishauser, A., Simonsen, S. B., Smith, A., Solgaard, A., Spencer, M., Steen-Larsen, H. C., Stevens, C. M., Sugiyama, S., Tedesco, M., Thompson-Munson, M., Tsutaki, S., Van_As, D., Van_Den_Broeke, M. R., Wilhelms, F., Xiao, J., and Xiao, C.: The SUMup collaborative database: Surface mass balance, subsurface temperature and density measurements from the Greenland and Antarctic ice sheets (1912 - 2023), https://doi.org/10.18739/A2M61BR5M, publisher: NSF Arctic Data Center, 2023.

970 Veldhuijsen, S. B. M., van de Berg, W. J., Brils, M., Kuipers Munneke, P., and van den Broeke, M. R.: Characteristics of the 1979–2020 Antarctic firn layer simulated with IMAU-FDM v1.2A, The Cryosphere, 17, 1675–1696, https://doi.org/10.5194/tc-17-1675-2023, 2023.

Villefranque, N., Blanco, S., Couvreux, F., Fournier, R., Gautrais, J., Hogan, R. J., Hourdin, F., Volodina, V., and Williamson, D.: Process-Based Climate Model Development Harnessing Machine Learning: III. The Representation of Cumulus Geometry and Their 3D Radiative Effects, Journal of Advances in Modeling Earth Systems, 13, e2020MS002423, https://doi.org/10.1029/2020MS002423, _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2020MS002423, 2021.

975 Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geoscientific Model Development, 5, 773–791, https://doi.org/10.5194/gmd-5-773-2012, publisher: Copernicus GmbH, 2012.

Walter, B., Weigel, H., Wahl, S., and Löwe, H.: Wind tunnel experiments to quantify the effect of aeolian snow transport on the surface snow microstructure, The Cryosphere, 18, 3633–3652, https://doi.org/10.5194/tc-18-3633-2024, publisher: Copernicus GmbH, 2024.

980 Wang, T., Ottlé, C., Boone, A., Ciais, P., Brun, E., Morin, S., Krinner, G., Piao, S., and Peng, S.: Evaluation of an improved intermediate complexity snow scheme in the ORCHIDEE land surface model, Journal of Geophysical Research: Atmospheres, 118, 6064–6079, https://doi.org/10.1002/jgrd.50395, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jgrd.50395, 2013.

Wang, Y., Ding, M., Wessem, J. M. v., Schlosser, E., Altnau, S., Broeke, M. R. v. d., Lenaerts, J. T. M., Thomas, E. R., Isaksson, E., Wang, J., and Sun, W.: A Comparison of Antarctic Ice Sheet Surface Mass Balance from Atmospheric Climate Models and In Situ Observations, Journal of Climate, 29, 5317–5337, https://doi.org/10.1175/JCLI-D-15-0642.1, publisher: American Meteorological Society Section: Journal of Climate, 2016.

Williamson, D., Goldstein, M., Allison, L., Blaker, A., Challenor, P., Jackson, L., and Yamazaki, K.: History matching for exploring and reducing climate model parameter space using observations and a large perturbed physics ensemble, Climate Dynamics, 41, 1703–1729, https://doi.org/10.1007/s00382-013-1896-4, 2013.
<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

990 Williamson, D. B., Blaker, A. T., and Sinha, B.: Tuning without over-tuning: parametric uncertainty quantification for the NEMO ocean model, Geoscientific Model Development, 10, 1789–1816, https://doi.org/10.5194/gmd-10-1789-2017, publisher: Copernicus GmbH, 2017.