Numerical experiments on vapor diffusion in polar snow and firn and its impact on isotopes using the multi-layer energy balance model Crocus in SURFEX v8.0

To evaluate the impact of vapor diffusion on isotopic composition variations in snow pits and then in ice cores, we introduced water isotopes in the detailed snowpack model Crocus. At each step and for each snow layer, (1) the initial isotopic composition of vapor is taken at equilibrium with the solid phase, (2) a kinetic fractionation is applied during transport, and (3) vapor is condensed or snow is sublimated to compensate for deviation to vapor pressure at saturation.

We study the different effects of temperature gradient, compaction, wind compaction, and precipitation on the final vertical isotopic profiles. We also run complete simulations of vapor diffusion along isotopic gradients and of vapor diffusion driven by temperature gradients at GRIP, Greenland and at Dome C, Antarctica over periods of 1 or 10 years. The vapor diffusion tends to smooth the original seasonal signal, with an attenuation of 7 to 12 % of the original signal over 10 years at GRIP. This is smaller than the observed attenuation in ice cores, indicating that the model attenuation due to diffusion is underestimated or that other processes, such as ventilation, influence attenuation. At Dome C, the attenuation is stronger (18 %), probably because of the lower accumulation and stronger δ¹⁸O gradients.

1 Introduction

The isotopic ratios of oxygen or deuterium measured in ice cores have been used for a long time to reconstruct the evolution of temperature over the Quaternary (EPICA comm. members, 2004; Johnsen et al., 1995; Jones et al., 2018; Jouzel et al., 2007; Kawamura et al., 2007; Lorius et al., 1985; Petit et al., 1999; Schneider et al., 2006; Stenni et al., 2004, 2011; Uemura et al., 2012; WAIS Divide project members, 2013). They are, however, subject to alteration during post-deposition through various processes. Consequently, even if the link between the temperature and isotopic composition of the precipitations is quantitatively determined from measurements and modeling studies (Stenni et al., 2016; Goursaud et al., 2017), it cannot faithfully be applied to reconstructions of past temperature. Nevertheless, ice cores remain a primary climatic archive for the Southern Hemisphere where continental archives are rare (Mann and Jones, 2003). In Antarctica, where meteorological records only started in the 1950s (Genthon et al., 2013), they provide useful information for understanding climate variability (e.g., EPICA comm. members, 2006; Shaheen et al., 2013; Steig, 2006; Stenni et al., 2011) and recent climate change (e.g., Altnau et al., 2015; Schneider et al., 2006). When using ice cores for past climate reconstruction, parameters other than temperature at condensation influence the isotopic compositions and must be considered. Humidity and temperature in the region of evaporation (Landais et al., 2008; Masson-Delmotte et al., 2011) or the seasonality of precipitation (Delmotte et al., 2000; Sime et al., 2008; Laepple et al., 2011) should be taken into account. In addition, uneven accumulation in time and space introduces stratigraphic noise (Ekaykin et al., 2009). Records from adjacent snow pits have been shown to be markedly different under the influence of decameter-scale local effects such as wind redeposition of snow, erosion, compaction, and metamorphism (Ekaykin et al., 2014; Petit et al., 1982). These local effects reduce the signal to noise ratio. Then only stacking a series of records from snow pits can eliminate this local variability and yield information relevant to recent climate variations (Fisher and Koerner, 1994; Hoshina et al., 2014; Ekaykin et al., 2014; Altnau et al., 2015). This concern is particularly significant in the central regions of East Antarctica characterized by accumulation rates lower than 100 mm of water equivalent per year (van de Berg et al., 2006). There, strong winds can scour and erode the snow layer over depths larger than the annual accumulation (Frezzotti et al., 2005; Morse et al., 1999; Libois et al., 2014). There is thus a strong need to study post-deposition effects in these cold and dry regions.

Additionally to the mechanical reworking of the snow, the isotopic compositions are further modified in the snowpack. First, diffusion along isotopic gradients can occur within the snow grains due to solid diffusion (Ramseier et al., 1967). Second, within the porosity, the vapor isotopic composition can change due to (1) diffusion along isotopic gradients in the gaseous state, (2) thermally induced vapor transport caused by vapor pressure gradients, (3) ventilation in the gaseous state, or (4) exchanges between the gas phase and the solid phase, i.e., sublimation and condensation. In the porosity, the combination of diffusion along isotopic gradients in the vapor and of exchange between vapor and the solid phase has been suggested to be the main explanation to the smoothing of the isotopic signal in the solid phase (Ebner et al., 2016, 2017; Gkinis et al., 2014; Johnsen et al., 2000). The isotopic compositions in the solid phase are also modified by “dry metamorphism” and “ventilation” but in a less predictable way. In both cases, the vapor transport exerts an influence on the isotopic compositions in the solid phase because of permanent exchanges between solid and vapor. During “dry metamorphism” (Colbeck et al., 1983), vapor transport is driven by vapor pressure gradients, themselves caused by temperature gradients. During ventilation (Town et al., 2008), vapor moves as part of the air in the porosity because of pressure variations at the surface. Last, at the top of the snowpack, the isotopic composition of snow may also be modified through direct exchange with atmospheric vapor (Ritter et al., 2016).

To elucidate the impact of these various post-deposition processes on the snow isotopic compositions, numerical models are powerful tools. They allow one to discriminate between processes and test their impact one at a time. For instance, Johnsen et al. (2000) were able to simulate and deconvolute the influence of diffusion along isotopic gradients in the vapor at two Greenland ice-core sites, GRIP and NGRIP, using a numerical model. To do this, they define a quantity called “diffusion length”, which is the mean displacement of a water molecule during its residence time in the porosity. Using a thinning model and an equation for the diffusivity of the water isotopes in snow, they compute this diffusion length as a function of depth. It is then used to compute the attenuation ratio A∕A_o and in the end retrieve the original amplitude A_o. Additionally, the effect of forced ventilation was investigated by Neumann (2003) and Town et al. (2008) using similar multi-layer numerical models. In these models, wind-driven ventilation forces atmospheric vapor into snow. There, the vapor is condensed, especially in layers colder than the atmosphere.

We focus on the movement of water isotopes in the vapor phase in the porosity, in the absence of macroscopic air movement. In that situation, the movement of vapor molecules in the porosity is caused by vapor pressure gradients or by diffusion along isotopic gradients. Note that in the first case, the vapor transport is “thermally induced”; i.e., the vapor pressure gradients directly result from temperature gradients within the snowpack. Thus, the first prerequisite for our model is to correctly simulate macroscopic energy transfer within the snowpack and energy exchange at the surface.

The transport of vapor molecules will affect the isotopic composition in the solid phase only if exchanges between vapor and solid are also implemented. Thus, the second prerequisite is that the model includes a description of the snow microstructure and of its evolution in time. Snow microstructure is typically represented by its emerging scalar properties such as density, specific surface area, and higher-order terms often referred to as “shape parameters” (e.g., Krol and Löwe, 2016). While the concept of “grain” bears ambiguity, it is a widely used term in snow science and glaciology that we employ here as a surrogate for “elementary microstructure element” without explicit reference to a formal definition, whether crystallographic or geometrical.

Crocus is a unidimensional multi-layer model of snowpack with a typically centimetric resolution initially dedicated to the numerical simulation of snow in temperate regions (Brun et al., 1992). It describes the evolution of the snow microstructure driven by temperature and temperature gradients during dry snow metamorphism using semi-empirical variables and laws. It has been used for ice-sheet conditions in polar regions, both Greenland and Antarctica (Brun et al., 2011; Lefebre et al., 2003; Fréville et al., 2013; Libois et al., 2014, 2015). In these regions, it gives realistic predictions of density and snow type profiles (Brun et al., 1992; Vionnet et al., 2012), snow temperature profile (Brun et al., 2011), and snow specific surface area and permeability (Carmagnola et al., 2014; Domine et al., 2013). It has been recently optimized for application to conditions prevailing at Dome C, Antarctica (Libois et al., 2014). This was necessary to account for specific conditions such as high snow density values at the surface and low precipitation amounts.

The Crocus model has high vertical spatial resolution and also includes an interactive simulation of snow metamorphism in near-surface snow and firn. Therefore, it is a good basis for the study of post-deposition effects in low accumulation regions. For the purposes of this study, we thus implemented vapor transport resulting from temperature gradients and the water isotope dynamics into the Crocus model. This article presents this double implementation and a series of sensitivity tests. A perfect match with observations is not anticipated, in part because not all relevant processes are represented in the model. This study thus represents a first step towards better understanding the impact of diffusion driven by temperature gradients on the snow isotopic composition.

2 Physical basis

The isotopic composition of the snow can evolve after deposition due to several processes. Here, we first give a brief overview of such processes at the macroscopic level. Section 2.1 thus deals with the modification of the isotopic composition of a centimetric to decametric snow layer after exchanges with the other layers. Second, we consider the evolution of the isotopic composition at the microscopic level, i.e., at the level of the microstructure. The macroscopic change in the isotopic composition results from both large-scale and small-scale processes. For instance, dry metamorphism includes both vapor transport from one layer to another and vapor–ice grain exchange inside a layer.

2.1 Evolution of the composition of the snow layers at the macroscopic scale

Several studies address the evolution of the isotopic compositions in the snow column after deposition. Here we describe first the processes leading only to attenuation of the original amplitude (Sect. 2.1.1). Then we describe processes that lead to other types of signal modifications (Sect. 2.1.2). These modifications result from the transportation and accumulation of heavy or light isotopes in some layers without any link to the original isotopic signal. In some cases, the mean δ¹⁸O value of the snow deposited can also be modified.

2.1.1 Signal attenuation on a vertical profile: smoothing

In this section we consider the processes leading only to attenuation of the original amplitude of the δ¹⁸O signal by smoothing. We define the mean local pluriannual value as the average isotopic composition in the precipitation taken over 10 years. The smoothing processes, which act only on signal variability, do not modify this average value. Within the snow layers, the smoothing of isotopic compositions is caused by diffusion along isotopic gradients in the vapor phase and solid phase. The magnitude of smoothing depends on site temperature and on accumulation. Higher temperatures correspond to higher vapor concentrations and higher diffusivities in the vapor and solid phases. Oppositely, high accumulation rates ensure a greater separation between seasonal δ¹⁸O peaks (Ekaykin et al., 2009; Johnsen et al., 1977), thereby limiting the impact of diffusion. They also result in increased densification rates and therefore reduced diffusivities (Gkinis et al., 2014). Because sites with high accumulation rates also usually have higher temperatures, the resulting effect on diffusion is still unclear. These two competing effects should be thoroughly investigated, and Johnsen et al. (2000) display the damping amplitude of a periodic signal depending on wavelength and on diffusion length, strongly driven by temperature.

In Greenland, Johnsen et al. (1977) indicate that annual cycles generally disappear at depths shallower than 100 m for sites with accumulation lower than 200 kg m⁻² yr⁻¹. Diffusion along isotopic gradients exists throughout the entire snow–ice column. It occurs mainly in the vapor phase in the firn, especially in the upper layers with larger porosities. After pore closure, it takes place mostly in the solid phase at a much slower rate. Note that in the solid phase, all isotopes have the same diffusion coefficient.

2.1.2 Signal shift caused by processes leading to oriented vapor transport

We consider here the oriented movement of water molecules forced by external variables such as temperature or pressure. We use the term “oriented” here to describe an overall movement of water molecules that is different from their molecular agitation and externally forced. Three processes can contribute to oriented vapor transport and hence possible isotopic modification within the snowpack: diffusion, convection, and ventilation (Albert et al., 2002). Brun and Touvier (1987) have demonstrated that the convection of dry air within the snow occurs only in the case of very low snow density of the order of ∼ 100 kg m⁻³. These conditions are generally not encountered in Antarctic snow and therefore convection is not considered here. Bartelt et al. (2004) also indicate that energy transfer by advection is negligible compared to energy transfer by conduction in the first meters of the snowpack. The two other processes of ventilation and diffusion are respectively forced by variations in the surface pressure and surface temperature. In the first case, the interaction between wind and surface roughness is responsible for wind pumping, i.e., the renewal of the air of the porosity through macroscopic air movement (Albert et al., 2002; Colbeck, 1989; Neumann and Waddington, 2004). In the second case, air temperature diurnal or seasonal variations generate vertical temperature gradients within the snow (Albert and McGilvary, 1992; Colbeck, 1983). They result in vertical vapor pressure gradients responsible for vapor diffusion. These two processes are largely exclusive (Town et al., 2008) because strong ventilation homogenizes the air and vapor in the porosity and therefore prevents diffusion. Diffusion as a result of temperature gradients can coexist with ventilation only at very low air velocities (Calonne et al., 2015). It becomes the main process of vapor transport when air is stagnant in the porosity. During diffusion, lighter molecules move more quickly in the porosity, leading to a kinetic fractionation of the various isotopologues.

2.2 Evolution of the isotopic composition at the microscopic scale

2.2.1 Conceptual representation of snow microstructure as spherical grains

The term “snow grain” as used classically is an approximation. In reality, snow grains are very diverse in size, shape, and degree of metamorphism and may also be made of several snow crystals agglomerated. Moreover, they are often connected to each other, forming an ice matrix, or “snow microstructure”. However, several studies addressing snow metamorphism physical processes have relied on spherical ice elements to represent snow grains and snow microstructure (Legagneux and Domine, 2005; Flanner and Zender, 2006). Here, we consider the snow grains to be made of two concentric layers, one internal and one external, with different isotopic compositions. In terms of snow microstructure, this could correspond to the inner vs. outer regions of the snow microstructure.

The snow grain or microstructure is not necessarily homogeneous in terms of isotopic composition. On the one hand, the central part of the grain or of the microstructure is relatively insulated. This central part becomes even more insulated as the grain grows or as the structure gets coarser. On the other hand, outer layers are not necessarily formed at the same time as the central part or in the same environment (Lu and DePaolo, 2016). They are prone to subsequent sublimation–condensation of water molecules, implying that their composition varies more frequently than for the inner layers. Of course, only the bulk δ¹⁸O value of the snow grain can be measured by mass spectrometry. But considering the heterogeneity of the grain may be required to get a fine understanding of the processes. In the following, we propose splitting the ice grain compartment into two sub-compartments: grain surface and grain center. Thus, the grain surface isotopic composition evolves because of exchange with two compartments: (1) vapor in the porosity through sublimation–condensation and (2) the grain center through solid diffusion or grain center translation. The grain center composition evolves at the timescale of weeks or months, as opposed to the grain surface, where the composition changes at the timescale of the vapor diffusion, i.e., over minutes.

2.2.2 Solid diffusion within snow grains

The grain center isotopic composition may change either as a result of crystal growth and/or sublimation or as a result of solid diffusion within the grain. For solid diffusion, water molecules move in the crystal lattice through a vacancy mechanism, in a process of self-diffusion that has no particular direction and that is very slow. The diffusivity of water molecules in solid ice D_ice in m² s⁻¹ follows the Arrhenius law. Thus, it can be expressed as a function of ice temperature T (Gkinis et al., 2014; Johnsen et al., 2000; Ramseier, 1967) using Eq. (1):

$$\begin{array}{}\text{(1)}& D_{\text{ice}}=\mathrm{9.2}\times {\mathrm{10}}^{-\mathrm{4}}\times \exp \left({\displaystyle \frac{-\mathrm{7186}}{T}}\right),\end{array}$$

for which the symbols are listed in Table 1.

Thus at 230 K, the diffusivity is $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{17}}$ m² s⁻¹. Gay et al. (2002) indicate that in the first meter at Dome C, a typical snow grain has a radius of 0.1 mm. Across this typical snow grain, the characteristic time for diffusion is given by Eq. (2).

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{t}_{\text{sol}}={\displaystyle \frac{R_{\text{moy}}^{\mathrm{2}}}{D_{\text{ice}}}}=\mathrm{4.03}\times {\mathrm{10}}^{\mathrm{8}}\text{s},\text{or}\sim \mathrm{13}\text{years}\end{array}$$

Table 1 Definition of the symbols used.

Therefore, the solid diffusion within the grain is close to zero at the timescales considered in the model. For Dome C, if we use the average temperature T of 248 K for the summer months (December to January; Table 2), the characteristic time becomes 15 months. Thus, within a summer period, the snow grain is only partially refreshed through this process. At Summit the grain size is typically larger, from 0.2 to 0.25 mm in wind-blown and wind-packed snow, and from 0.5 to 2 mm in the depth hoar layer (Albert and Shultz, 2002). The summer temperature is also higher, with an average value T of 259 K at Summit from July to September, after Shuman et al. (2001). Using a grain size of 0.25 mm, the resulting characteristic time is of the order of 30 months.

2.2.3 Snow grain recrystallization

During snow metamorphism, the number of snow grains tends to decrease with time, while the snow grain size tends to increase (Colbeck, 1983). Each grain experiences continuous recycling through sublimation–condensation, but the small grains are more likely to disappear completely. Then, there is no more nucleus for condensation at the grain initial position. Oppositely, the bigger grains do not disappear and accumulate the vapor released by the smaller ones. Concurrently to this change in grain size, the grain shape also tends to evolve. In conditions of a maintained or stable temperature gradient, facets appear at the condensing end of snow grains, while the sublimating end becomes rounded (Colbeck, 1983). In that case, the center of the grain moves toward the warm air region. This migration causes a renewal of the grain center on a proportion that can be estimated from the apparent grain displacement (Pinzer et al., 2012). Pinzer et al. (2012) use this method to obtain an estimation of vapor fluxes.

The asymmetric recrystallization of snow grains implies that the surface layer of the snow grain is eroded at one end and buried at the other end. Therefore, the composition of the grain center changes more often than if the surface layer was thickening through condensation or thinning through sublimation homogeneously over the grain surface. This means that the “inner core” of the grain gets exposed more often. Implementing this process is thus very important to have a real-time evolution of the snow grain center isotopic composition. Here, we reverse the method of Pinzer et al. (2012). Therefore, we use the fluxes of isotopes in the vapor phase computed by the model to assess the renewal of the grain center (Sect. 3.1.3.).

3 Material and methods

3.1 Description of the model SURFEX–Crocus v8.0

We first present the model structure and second describe the new module of vapor transport (diffusion forced by temperature gradients). Third, we present the integration of water isotopes in the model.

3.1.1 Model structure

The Crocus model is a one-dimensional detailed snowpack model consisting of a series of snow layers with variable and evolving thicknesses. Each layer is characterized by its density, heat content, and by parameters describing snow microstructure such as sphericity and specific surface area (Vionnet et al., 2012; Carmagnola et al., 2014). In the model, the profile of temperature evolves with time in response to (1) surface temperature and (2) energy fluxes at the surface and at the base of the snowpack. To correctly compute energy balance, the model integrates albedo calculation as deduced from surface microstructure and impurity content (Brun et al., 1992; Vionnet et al., 2012).

The successive components of the Crocus model have been described by Vionnet et al. (2012). Here we only list them to describe those modified to include water stable isotopes and water vapor transfer. Note that the Crocus model has a typical internal time step of 900 s (15 min) corresponding to the update frequency of layers properties. We only refer here to processes occurring in dry snow.

1. Snowfall. The presence or absence of precipitation at a given time is determined from the atmospheric forcing inputs. When there is precipitation, a new layer of snow may be formed. Its thickness is deduced from the precipitation amount.

2. Update of snow layering. At each step, the model may split one layer into two or merge two layers together to get closer to a target vertical profile for optimal calculations. This target profile has high resolution in the first layers to correctly simulate heat and matter exchanges. The layers that are merged together are the closest in terms of microstructure variables.

3. Metamorphism. The evolution of microstructure variables follows empirical laws. These laws describe the change in grain parameters as a function of temperature, temperature gradient, snow density, and liquid water content.

4. Snow compaction. Layer thickness decreases and layer density increases under the burden of the overlying layers and resulting from metamorphism. In the original module, snow viscosity is parameterized using the layer density and also using information on the presence of hoar or liquid water. However, this parameterization of the viscosity was designed for alpine snowpack (Vionnet et al., 2012) and may not be adapted to polar snowpacks. Moreover, since we are considering only the first 12 m of the snowpack in the present simulations, the compaction in the considered layers does not compensate for the yearly accumulation, leading to a rising snow level with time. To maintain a stable surface level in our simulations, we used a simplified compaction scheme in which the compaction rate ε is the same for all the layers. The compaction rate is obtained by dividing the accumulation rate at the site (see Sect. 3.3) by the total mass of the snow column (Eq. 3). It is then applied to all layers to obtain the density change per time step using Eq. (4).

   $$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{\epsilon }={\displaystyle \frac{\text{d}{m}_{\text{sn}}}{\text{d}t}}/\sum _{\mathrm{1}}^{\text{nmax}}\left({\mathit{\rho }}_{\text{sn}}\left(t,n\right)\times \text{d}z(t,n)\right)\text{(4)}& {\displaystyle }& {\displaystyle \frac{{\mathit{\rho }}_{\text{sn}}(t+\text{d}t,n)-{\mathit{\rho }}_{\text{sn}}(t,n)}{\text{d}t}}=\mathit{\epsilon }\times {\mathit{\rho }}_{\text{sn}}(t,n)\end{array}$$

5. Wind drift events. They modify the properties of the snow grains, which tend to become more rounded. They also increase the density of the first layers through compaction. An option allows snow to be partially sublimated during these wind drift events (Vionnet et al., 2012).

6. Snow albedo and transmission of solar radiation. In the first 3 cm of snow, snow albedo and absorption coefficient are computed from snow microstructure properties and impurity content. The average albedo value in the first 3 cm is used to determine the part of incoming solar radiation reflected at the surface. The rest of the radiation penetrates into the snowpack. Then, the absorption coefficient is used to describe the rate of decay of the radiation as it is progressively absorbed by the layers downward, following an exponential law.

7. Latent and sensible surface energy and mass fluxes. The sensible heat flux and the latent heat flux are computed using the aerodynamic resistance and the turbulent exchange coefficients.

8. Vertical snow temperature profile. It is deduced from the heat diffusion equation using the snow conductivity and the energy balance at the top and at the bottom of the snowpack.

9. Snow sublimation and condensation at the surface. The amount of snow sublimated–condensed is deduced from the latent heat flux, and the thickness of the first layer is updated. Other properties of the first layer such as density and SSA are kept constant.

3.1.2 Implementation of water transfer

The new vapor transport subroutine has been inserted after the compaction (Eq. 4) and wind drift (Eq. 5) modules and before the solar radiation module (Eq. 6). In this section, the term “interface” is used for the horizontal surface of exchange between two consecutive layers. The flux of vapor at the interface between two layers is obtained using Fick's law of diffusion (Eq. 5):

$$\begin{array}{ll}\text{(5)}& {\displaystyle }& {\displaystyle }F(n+\mathrm{1}\to n)={\displaystyle }& {\displaystyle \frac{-\mathrm{2}{D}_{\text{eff}}\left(t,n\to n+\mathrm{1}\right)\left(C_{\text{v}}(t,n)-C_{\text{v}}(t,n+\mathrm{1})\right)}{\text{d}z(t,n)+\text{d}z(t,n+\mathrm{1})}},\end{array}$$

where dz(t,n) and $\text{d}z(t,n+\mathrm{1})$ are the thicknesses of the two layers considered in meters, C_v(t,n) and $C_{\text{v}}(t,n+\mathrm{1})$ are the local vapor mass concentrations in the two layers in kg m⁻³, and $D_{\text{eff}}\left(t,n\to n+\mathrm{1}\right)$ in m² s⁻¹ is the effective diffusivity of water vapor in the snow at the interface. The thicknesses are known from the previous steps of the Crocus model, but the vapor mass concentrations and the interfacial diffusivities must be computed.

The effective diffusivity at the interface is obtained in two steps: first the effective diffusivities (D_eff(t,n) and $D_{\text{eff}}(t,n+\mathrm{1})$) in each layer are calculated (Eq. 6); second, the interfacial diffusivity ($D_{\text{eff}}(t,n\to n+\mathrm{1}))$ is computed as their harmonic mean (Eq. 7). Effective diffusivity can be expressed as a function of the snow density using the relationship proposed by Calonne et al. (2014) for layers with relatively low density. In these circumstances, the compaction occurs by “boundary sliding”, meaning that the grains slide on each other, but that their shape is not modified. It is therefore applicable to our study for which density is always below 600 kg m⁻³. The equation of Calonne et al. (2014) is based on the numerical analysis of 3-D tomographic images of different types of snow. It relates normalized effective diffusivity D_eff∕D_v to the snow density ρ_sn in the layer (Eq. 6). D_v is the vapor diffusivity in air and has a value that varies depending on the air pressure and air temperature (Eq. 19 in Johnsen et al., 2000). ρ_ice corresponds to the density of ice and has a value of 917 kg m⁻³.

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle \frac{D_{\text{eff}}(t,n)}{D_{\text{v}}}}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\text{sn}}(t,n)}{{\mathit{\rho }}_{\text{ice}}}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\text{(7)}& {\displaystyle }& {\displaystyle }{D}_{\text{eff}}\left(t,n\to n+\mathrm{1}\right)={\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{1}}{D_{\text{eff}}\left(t,n\right)}+\frac{\mathrm{1}}{D_{\text{eff}}\left(t,n+\mathrm{1}\right)}}}\end{array}$$

We assume that vapor is generally at saturation in the snow layers (Neumann et al., 2008, 2009). The local mass concentration of vapor C_v in kg m⁻³ in each layer is given by the Clausius–Clapeyron equation (Eq. 8):

$$\begin{array}{}\text{(8)}& C_{\text{v}}(t,n)=C_{\text{v0}}\exp \left({\displaystyle \frac{L_{\text{sub}}}{R_{\text{v}}{\mathit{\rho }}_{\text{ice}}}}\left({\displaystyle \frac{\mathrm{1}}{T_{\mathrm{0}}}}-{\displaystyle \frac{\mathrm{1}}{T(t,n)}}\right)\right),\end{array}$$

where C_v0 is the mass concentration of vapor at 273.16 K and is equal to $\mathrm{2.173}\times {\mathrm{10}}^{-\mathrm{3}}$ kg m⁻³, L_sub is the latent heat of sublimation and has a value of 2.6×10⁹ J m⁻³, R_v is the vapor constant and has a value of 462 J kg⁻¹ K⁻¹, ρ_ice is the density of ice and has a value of 917 kg m⁻³, T₀ is the temperature of the triple point of water and is equal to 273.16 K, and T is the temperature of the layer.

All layers are treated identically, except the first layer at the top and the last layer at the bottom. For the uppermost layer, the exchange of vapor occurs only at the bottom boundary. Exchanges with the atmosphere are described elsewhere in Crocus at step 9 in which surface energy balance is realized. For the lowermost layer, only exchanges taking place at the top boundary are considered, with the flux of vapor to and from the underlying medium being set to zero.

For each layer, the mass concentration of vapor in the air and the effective diffusivity are computed within the layer and in the neighboring layers. Fluxes at the top and bottom of each layer are deduced from Fick's law of diffusion (Eq. 5). They are integrated over the subroutine time step, and the new mass of the layer is computed. It is used at the beginning of the next subroutine step. We use a 1 s time step within the subroutine, which is smaller than the main routine time step of 900 s. This ensures that vapor fluxes remain small relative to the amount of vapor present in the layers. Note that the temperature profile, which controls the vapor pressure profile, is not modified within the subroutine. Physically, temperature values should change as a result of the transfer of sensible heat from one layer to another associated with vapor transport. They should also evolve due to the loss or gain of heat caused by water sublimation–condensation (Albert and McGilvary, 1992; Kaempfer et al., 2005). However, vapor transport is only a small component of heat transfer between layers (Albert and Hardy, 1995; Albert and McGilvary, 1992). In the absence of ventilation, with or without vapor diffusion, the steady-state profile for temperature varies by less than 2 % (Calonne et al., 2014). Thus, the effect can be neglected at first order.

3.1.3 Implementation of water isotopes

In the model, the isotopic composition of snow in each layer is represented by the triplicate (δ¹⁸O, d-excess, ¹⁷O-excess). Only the results of δ¹⁸O are presented and discussed here. For each parameter, two values per layer are considered independently that correspond to the “snow grain center” and the “snow grain surface”. Water vapor isotopic composition is deduced at each step from the snow grain surface isotopic composition. It is not stored independently to limit the number of prognostic variables. The isotopic compositions are used at step 1, i.e., for snowfall, and after step 5 within the new module of vapor transfer.

In the snowfall subroutine, a new layer of snow may be added, depending on the weather, at the top of the snowpack. At this step of the routine, the snow grains being deposited are presumed to be homogenous, i.e., they have the same composition in the grain surface compartment and in the grain center compartment. Their composition is deduced from the air temperature (see Sect. 3.2).

Within the vapor transport subroutine, a specific module deals with the isotopic aspects of vapor transport. It modifies the isotopic compositions in the two snow grain sub-compartments as a result of water vapor transport and the recrystallization of snow crystals. It works with four main steps:

1. an initiation step in which the vapor isotopic compositions are computed using equilibrium fractionation from the ones in the grain surface sub-compartment;

2. a transport step in which vapor moves from one layer to another, with a kinetic fractionation associated with diffusion;

3. a balance step in which the new vapor in the porosity exchanges with the grain surface compartment by sublimation–condensation (the flux is determined by the difference between the actual vapor mass concentration and the expected vapor mass concentration at saturation); and

4. a “recrystallization” step in which the grain center and grain surface isotopic compositions are homogenized, leading to an evolution of grain center isotopic composition.

The time step in this module is 1 s, the same as the time step of the subroutine.

The initial vapor isotope composition $R_{\text{vap,ini}}^i$ in a given layer is taken at equilibrium with the grain surface isotopic composition $R_{\text{surf,ini}}^i$. Here i denotes a heavy isotope and thus stands for ¹⁸O, ¹⁷O, or D. Equilibrium fractionation is a hypothesis that is correct in layers in which vapor has reached equilibrium with ice grains both physically and chemically. This process is limited by the water vapor–snow mass transfer whose associated speed is of the order of 0.09 m s⁻¹ (Albert and McGilvary, 1992). In our case, we are dealing with centimetric-scale layer thickness and recalculate the isotopic composition every second so that we consider the speed of the mass transfer as not limiting the equilibrium situation at the water vapor–snow interface. To compute isotopic ratios for water vapor we use the following Eqs. (9) and (10).

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{c}{R}_{\text{vap,ini}}^{\mathrm{18}}={\mathit{\alpha }}_{\text{sub}}^{\mathrm{18}}\times R_{\text{surf,ini}}^{\mathrm{18}}\\ R_{\text{vap,ini}}^{\mathrm{17}}={\mathit{\alpha }}_{\text{sub}}^{\mathrm{17}}\times R_{\text{surf,ini}}^{\mathrm{17}}\\ R_{\text{vap,ini}}^{\mathrm{18}}+R_{\text{vap,ini}}^{\mathrm{17}}+\mathrm{1}=\mathrm{1}/c_{\text{vap,ini}}^{\mathrm{16}}\end{array}\right.\text{(10)}& {\displaystyle }& {\displaystyle }\left\{\begin{array}{c}{R}_{\text{vap,ini}}^D={\mathit{\alpha }}_{\text{sub}}^D\times R_{\text{surf,ini}}^D\\ R_{\text{vap,ini}}^D+\mathrm{1}=\mathrm{1}/c_{\text{vap,ini}}^{\mathrm{1}H}\end{array}\right.\end{array}$$

The equilibrium fractionation coefficients (${\mathit{\alpha }}_{\text{sub}}^i)$ are obtained using the temperature-based parameterization from Ellehoj et al. (2013). Note that we make a slight approximation here by replacing molar concentrations with mass concentrations in our mass balance formulas (see Table 1 for symbol definitions).

The initial vapor mass concentration in air C_v has already been computed in the vapor transport subroutine, and the volume of the porosity can be obtained from the snow density ρ_sn and the thickness of the layer dz. By combining both, we obtain Eq. (11), which gives the initial mass of vapor in the layer m_vap,ini.

$$\begin{array}{}\text{(11)}& m_{\text{vap,ini}}=C_{\text{v}}\times \left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\text{sn}}}{{\mathit{\rho }}_{\text{ice}}}}\right)\times \text{d}z\end{array}$$

This mass of vapor should be subtracted from the initial grain surface mass because vapor mass is not tracked outside of the subroutine (Fig. 1). The new grain surface isotope composition, after vapor individualization is given by Eq. (12).

$$\begin{array}{}\text{(12)}& c_{\text{surf,new}}^{\mathrm{18}}={\displaystyle \frac{m_{\text{surf,new}}^{\mathrm{18}}}{m_{\text{surf,new}}}}={\displaystyle \frac{m_{\text{surf,ini}}^{\mathrm{18}}-m_{\text{vap,ini}}\times c_{\text{vap,ini}}^{\mathrm{18}}}{m_{\text{surf,ini}}-m_{\text{vap,ini}}}}\end{array}$$

The diffusion of isotopes follows the same scheme as the water vapor diffusion described above in Sect. 3.1.2. and Eq. (5). In Eq. (13), the gradient of vapor mass concentrations is replaced by a gradient of concentration for the studied isotopologue. The kinetic fractionation during the diffusion is realized with the Dⁱ∕D term where i stands for ¹⁸O, ¹⁷O, or ²H (Barkan and Luz, 2007).

$$\begin{array}{ll}{\displaystyle }{F}^{\mathrm{18}}& {\displaystyle }(n+\mathrm{1}\to n)={\displaystyle \frac{-\mathrm{2}\times D_{\text{eff}}\left(t,n\to n+\mathrm{1}\right)}{\text{d}z(t,n)+\text{d}z(t,n+\mathrm{1})}}\times {\displaystyle \frac{D^{\mathrm{18}}}{D}}\\ {\displaystyle }& {\displaystyle }\times \left(C_{\text{v}}(t,n)\times c_{\text{vap,ini}}^{\mathrm{18}}(t,n)-C_{\text{v}}(t,n+\mathrm{1})\right.\\ \text{(13)}& {\displaystyle }& {\displaystyle }\left.\times c_{\text{vap,ini}}^{\mathrm{18}}(t,n+\mathrm{1})\right)\end{array}$$

As done for water molecule transport (Sect. 3.1.2), the flux is set to zero at the top of the first layer and at the bottom of the last layer. When the vapor concentration is the same in two adjacent layers, the total flux of vapor is null. But diffusion along isotopic gradients still occurs if the isotopic gradients are nonzero (Eq. 13). Once the top and bottom fluxes of each layer have been computed, the new masses of the various isotopes in the vapor are deduced, as are the new ratios.

Figure 1 Splitting of the snow layer into two compartments, grain center and grain surface, with a constant mass ratio between them. The vapor compartment is a sub-compartment inside the grain surface compartment and is only defined at specific steps of the model.

After the exchanges between layers, the isotopic composition in the vapor has changed. However, the vapor isotopic composition is not a prognostic variable outside of the vapor transport subroutine. To record this change, it must be transferred to either the grain surface compartment or to the grain center compartment before leaving the subroutine. First, we consider exchanges of isotopes with the grain surface compartment, which is in direct contact with the vapor. Depending on the net mass balance of the layer, two situations must be considered.

1. If the mass balance is positive, condensation occurs so that the transfer of isotopes takes place from the vapor toward the grain surface. To evaluate the change in the isotope composition in the grain surface, the mass of vapor condensed Δm_vap,exc must be computed. It is the difference between the mass of vapor expected at saturation and the mass of vapor present in the porosity after vapor transport. Note that temperature does not evolve in this subroutine. Nevertheless, the difference is not exactly equal to the mass of vapor that has entered the layer because of layer porosity change. The excess mass of vapor is given by Eq. (14).

   $$\begin{array}{ll}\text{(14)}& {\displaystyle }\mathrm{\Delta }{m}_{\text{vap,exc}}& {\displaystyle }=\left[\left({\mathit{\rho }}_{\text{sn,new}}-{\mathit{\rho }}_{\text{sn,ini}}\right)+C_{\text{v}}\right.{\displaystyle }& {\displaystyle }\left.\times \left[\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\text{sn,ini}}}{{\mathit{\rho }}_{\text{ice}}}}\right)-\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\text{sn,new}}}{{\mathit{\rho }}_{\text{ice}}}}\right)\right]\right]\times \text{d}z\end{array}$$

   Since the excess of vapor is positive, the next step is the condensation of the excess vapor. The number of excess water molecules is determined through comparison with the expected number in the water vapor phase for an equilibrium state between surface snow and water vapor. Here the condensation of excess vapor occurs without additional fractionation because (1) there is a permanent isotopic equilibrium between surface snow and interstitial vapor restored at each first step of the subroutine, and (2) kinetic fractionation associated with diffusion is taken into account during the diffusion of the different isotopic species along the isotopic gradients.

2. If the mass balance is negative, the transfer of isotopes takes place from the grain surface toward the vapor without fractionation. Ice from the grain surface sub-compartment is sublimated without fractionation to reach the expected vapor concentration at saturation. Note that the absence of fractionation at sublimation is a frequent hypothesis because water molecules move very slowly in ice lattice (Friedman et al., 1991; Neumann et al., 2005; Ramseier, 1967). Consequently, the sublimation removes all the water molecules present at the surface of grains, including the heaviest ones, before accessing inner levels. In reality, there is evidence for fractionation at sublimation. It occurs through kinetic effects associated with sublimation or simultaneous condensation or during equilibrium fractionation at the boundary, especially when invoking the existence of a thin liquid layer at the snow–air interface (Neumann et al., 2008 and references therein; Sokratov and Golubev, 2009; Stichler et al., 2001; Ritter et al., 2016). The new composition in the vapor results from a mixing between the vapor present and the new vapor recently produced. The composition in the grain surface ice compartment does not change.

The limit between the surface compartment and the grain center compartment is defined by the mass ratio of the grain surface compartment to the total grain mass, i.e., $\mathit{\tau }=m_{\text{surf}}/(m_{\text{center}}+m_{\text{surf}}$); see also Fig. 1. This mass ratio can be used to determine the thickness of the grain surface layer as a fraction of grain radius for spherical grains. The surface compartment must be thin to be able to react to very small changes in mass when vapor is sublimated–condensed. Our model has a numerical precision of 6 decimals and is run at a 1 s temporal resolution. Consequently, the isotopic composition of the surface compartment can change in response to surface fluxes only if its mass is smaller than 10⁶ times the mass of the water vapor present in the porosity. This constrains the maximum value for τ: m_surf<10⁶m_vap or $m_{\text{surf}}/(m_{\text{center}}+m_{\text{surf}})<{\mathrm{10}}^{\mathrm{6}}\frac{\mathrm{\Phi }\cdot {\mathit{\rho }}_{\text{v}}\cdot V_{\text{tot}}}{{\mathit{\rho }}_{\text{sn}}\cdot V_{\text{tot}}}$, i.e., $\mathit{\tau }<\frac{{\mathit{\rho }}_v\cdot \mathrm{\Phi }}{{\mathit{\rho }}_{\text{sn}}}\times {\mathrm{10}}^{\mathrm{6}}$. Considering typical temperatures, snow densities, and layer thicknesses (Table 3), we obtain a maximum value of 3.3 $\times {\mathrm{10}}^{-\mathrm{2}}$. On the other hand, this compartment must be thick enough to transmit the change in isotopic compositions caused by vapor transport and condensation–sublimation to the grain center. Again, numerical precision imposes that its mass should be no less than 10⁻⁶ times the mass of the grain center compartment, and thus we get an additional constraint: $\mathit{\tau }>{\mathrm{10}}^{-\mathrm{6}}$. Here we use a ratio $\mathit{\tau }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ for the mass of the grain surface relative to the total mass of the layer (Fig. 1). We have run sensitivity tests with smaller and larger ratios (Sect. 4.3).

Two types of mixing between the grain surface and grain center are implemented in the model. The first one is associated with crystal growth or shrinkage because of vapor transfer. Mixing is performed at the end of the vapor transfer subroutine after sublimation–condensation has occurred. During the exchange of water between vapor and the grain surface, the excess or default of mass in the water vapor caused by vapor transport has been entirely transferred to the grain surface sub-compartment. Thus, the mass ratio between the grain surface compartment and the grain center compartment deviates from the original one. To bring the ratio τ back to the normal value of 5 $\times {\mathrm{10}}^{-\mathrm{4}}$, mass is transferred either from the grain surface to the grain center or from the grain center to the grain surface. This happens without fractionation; i.e., if the transfer occurs from the center to the surface, the composition of the center remains constant.

The second type of mixing implemented is the grain center translation (Pinzer et al., 2012), which favors mixing between the grain center and grain surface in the case of a sustained temperature gradient. Pinzer et al. (2012) used the apparent grain displacement to compute vapor fluxes. Here, we reverse this method and use the vapor fluxes computed from Fick's law to estimate the grain center renewal. We could transfer a small proportion of the surface compartment to the grain center every second. Instead, we choose to totally mix the snow grain every few days. The interval Δt_surf/center between two successive mixings is derived from the vapor flux $F(n+\mathrm{1}\to n)$ within the layer using Eq. (15).

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }{t}_{\text{surf/center}}={\displaystyle \frac{m_{\text{sn}}\times \mathit{\tau }}{F(n+\mathrm{1}\to n)}}\end{array}$$

The average temperature gradient of 3 ^∘C m⁻¹ corresponds to a flux $F(n+\mathrm{1}\to n)$ of $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{9}}$ kg m⁻² s⁻¹. The typical mass for the layer m_sn is 3.3 kg. Based on these values, the dilution of the grain surface compartment into the grain center should occur every 15 days. Of course, this is only an average, since layers have varying masses and since the temperature gradient can be larger or smaller. We will, however, apply this time constant for all the layers and any temperature gradient (see sensitivity tests Sect. 4.3) to ensure that the mixing between compartments occurs at the same time in all layers.

In terms of magnitude, this process is probably much more efficient for mixing the solid grain than grain growth or solid diffusion. It is thus crucial for the modification of the bulk isotopic composition of the snow layer. It creates the link between microscopic processes and macroscopic results.

Table 2 Climate and isotope variability at GRIP (Greenland) and Dome C (Antarctica).

Table 3 Typical thickness, density, temperature, and other parameters of the snow layers in the simulations. The ratio τ is the mass ratio between the grain surface compartment and the grain center compartment. It must be chosen within the interval [10⁻⁶; 10⁶ (C_vΦ∕ρ_sn)] to allow for exchanges between the grain surface compartment and grain center compartment on the one hand and between the grain surface compartment and vapor compartment on the other hand (see text for details).

3.1.4 Model initialization

For model initialization, an initial snowpack is defined with a fixed number of snow layers and for each snow layer an initial value of thickness, density, temperature, and δ¹⁸O. Typically, processes of oriented vapor transport such as thermally induced diffusion and ventilation occur mainly in the first meters of snow. Therefore, the model starts with an initial snowpack of about 12 m.

The choice of the layer thicknesses depends on the annual accumulation. Because the accumulation is much higher at GRIP than at Dome C (Sect. 3.2., Table 2), the second site is used to define the layer thicknesses. About 10 cm of fresh snow is deposited every year (Genthon et al., 2016; Landais et al., 2017). This implies that to keep seasonal information, at least one point every 4 cm is required in the first meter. For the initial profile, we impose a maximal thickness of 2 cm for the layers between 0 and 70 cm of depth and 4 cm for the layers between 70 cm and 2 m of depth. As the simulation runs, merging is allowed but restricted in the first meter to a maximum thickness of 2.5 cm. Below 2 m, the thicknesses are set to 40 cm or even 80 cm. Thus, the diffusion process can only be studied in the first 2 m of the model snowpack. In the very first centimeters of the snowpack, thin millimetric layers are used to accommodate low precipitation amounts and surface energy balance. The initial density profiles are defined for each site specifically (see Sect. 3.2). The initial temperature and δ¹⁸O profiles in the snowpack depend on the simulation considered (see Sect. 3.3).

3.1.5 Model output

A data file containing the spatiotemporal evolution of prognostic variables such as temperature, density, SSA, and δ¹⁸O is produced for each simulation. Here, we present the results for each variable as two-dimensional graphs, with time on the horizontal axis and snow height on the vertical axis. The variations in the considered variable are displayed as color levels. The white color corresponds to an absence of change in the variable. As indicated above, only the first 12 m of the polar snowpack are included in the model. The bottom of this initial snowpack constitutes the vertical reference or “zero” to measure vertical heights h. The height of the top of the snowpack varies with time due to snow accumulation and snow compaction. In the text, we sometimes refer to the layer depth z instead of its height h. The depth can be computed at any time by subtracting the current height of the considered layer from the current height of the top of the snowpack.

3.2 Studied sites: meteorology and snowpack description

In this study we run the model under the conditions encountered at Dome C, Antarctica and GRIP, Greenland. We chose these two sites because they have been well studied in recent years through field campaigns and numerical experiments. In particular for Dome C, a large amount of meteorological and isotopic data is available (Casado et al., 2016a; Stenni et al., 2016; Touzeau et al., 2016). Typical values of the main climatic parameters for the two studied sites, GRIP and Dome C, are given in Table 2 along with the typical δ¹⁸O range. Dome C has lower accumulation rates of 2.7 cm ice equivalent per year (ice eq. yr⁻¹) compared to GRIP rates of 23 cm ice eq. yr⁻¹ (Table 2), making it more susceptible to post-deposition processes.

Table 4 List of simulations described in the article with the corresponding paragraph number. The external atmospheric forcing used for Dome C is ERA-Interim reanalysis (2000–2013). However, the precipitation amounts from the ERA-Interim reanalysis are increased by 1.5 times to account for the dry bias in the reanalysis (as in Libois et al., 2014). For the second simulation at GRIP, Greenland meteorological conditions are derived from the atmospheric forcing of Dome C, but the temperature is modified ($T_{\text{GRIP}}=T_{\text{DC}}+\mathrm{15}$) as is the longwave down (LW_GRIP=0.85 LW_DC+60).

¹ Using data from 1-year snowfall sampling at Dome C (Stenni et al., 2016; Touzeau et al., 2016), we obtained the following Eq. (16) linking δ¹⁸O in the snowfall to the local temperature: ${\mathit{\delta }}^{\mathrm{18}}{O}_{\text{sf}}=\mathrm{0.45}\times (T-\mathrm{273.15})-\mathrm{31.5}$.
² The exponential profile of temperature used in Simulation 6 is defined using Eq. (20): $T\left(z\right)=T\left(\text{10\hspace{0.17em}m}\right)+\mathrm{\Delta }T\times \exp \left(-z/z\mathrm{0}\right)+\mathrm{0.1}\times z$
with T(10 m)=218 K, ΔT=28 K, and z0=1.516 m. It fits well with temperature measurements of midday in January (Casado et al., 2016b).
³ The Greenland snowpack has an initial sinusoidal profile of δ¹⁸O defined using Eq. (19): ${\mathit{\delta }}^{\mathrm{18}}O=-\mathrm{35.5}-\mathrm{8}\times \sin \left(\frac{\mathrm{2}\mathit{\pi }\times z}{a_{\text{GR}}\times {\mathit{\rho }}_{\text{ice}}/{\mathit{\rho }}_{\text{sn}}}\right)$.

Table 5 List of the sensitivity tests performed at GRIP and at Dome C. The external atmospheric forcing used for Dome C is ERA-Interim reanalysis (see Table 4).

¹The Greenland snowpack has an initial sinusoidal profile of δ¹⁸O defined using Eq. (19): ${\mathit{\delta }}^{\mathrm{18}}\text{O}=-\mathrm{35.5}-\mathrm{8}\times \sin \left(\frac{\mathrm{2}\mathit{\pi }\times z}{a_{\text{GR}}\times {\mathit{\rho }}_{\text{ice}}/{\mathit{\rho }}_{\text{sn}}}\right).$
² The Dome C snowpack has an initial sinusoidal profile of δ¹⁸O defined using Eq. (21): ${\mathit{\delta }}^{\mathrm{18}}\text{O}=-\mathrm{48.5}-\mathrm{6.5}\times \sin \left(\frac{\mathrm{2}\mathit{\pi }\times z}{a_{\text{DC}}\times {\mathit{\rho }}_{\text{ice}}/{\mathit{\rho }}_{\text{sn}}}\right).$

In this study, we also compare the results obtained for GRIP to results from two other Greenland sites, namely NGRIP and NEEM. GRIP is located at the ice-sheet summit, whereas the two other sites are located further north in lower elevation areas with higher accumulation rates. NGRIP is located 316 km to the NNW of the GRIP ice-drill site (Dahl-Jensen et al., 1997). GRIP and NGRIP have similar temperatures of −31.6 and −31.5 ^∘C but different accumulation rates of 23 and 19.5 cm ice eq. yr⁻¹, respectively. The NEEM ice-core site is located some 365 km to the NNW of NGRIP on the same ice ridge. It has an average temperature of −22 ^∘C and an accumulation rate of 22 cm ice eq. yr⁻¹.

The δ¹⁸O value in the precipitation at a given site reflects the entire history of the air mass, including evaporation, transport, distillation, and possible changes in trajectory and sources. However, assuming that these processes are more or less repeatable from one year to the next, it is possible to empirically relate the δ¹⁸O to the local temperature using measurements from collected samples. Here, using data from 1-year snowfall sampling at Dome C (Stenni et al., 2016; Touzeau et al., 2016), we use the following Eq. (16) to link δ¹⁸O in the snowfall to the local temperature T_air in K.

$$\begin{array}{}\text{(16)}& {\mathit{\delta }}^{\mathrm{18}}{\text{O}}_{\text{sf}}=\mathrm{0.45}\times \left(T_{\text{air}}-\mathrm{273.15}\right)-\mathrm{31.5}\end{array}$$

We do not provide an equivalent expression for GRIP, Greenland because the simulations run here (see Sect. 3.1.1) do not include precipitation.

The initial density profile in the snowpack is obtained from fitting density measurements from Greenland and Antarctica (Bréant et al., 2017). Over the first 12 m of snow, we obtain the following evolution (Eqs. 17 and 18) for GRIP and Dome C, respectively.

$$\begin{array}{ll}\text{(17)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\text{sn}}(t,n)=\mathrm{17.2}\times z(t=\mathrm{0},n)+\mathrm{310.3}{\displaystyle }& {\displaystyle }(N=\mathrm{22};R^{\mathrm{2}}=\mathrm{0.95})\\ \text{(18)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\text{sn}}(t,n)=\mathrm{12.41}\times z(t=\mathrm{0},n)+\mathrm{311.28}{\displaystyle }& {\displaystyle }(N=\mathrm{293};R^{\mathrm{2}}=\mathrm{0.50})\end{array}$$

3.3 List of simulations

Table 4 presents the model configuration for the six simulations considered here. Additionally, Table 5 presents sensitivity tests performed to evaluate the uncertainties associated with grain renewal parameters.

3.3.1 Greenland simulations

The first simulation, listed as number 1 in Table 4, is dedicated to the study of diffusion along isotopic gradients. It is realized on a Greenland snowpack with an initial sinusoidal profile of δ¹⁸O (see Eq. 19) and with a uniform and constant vertical temperature profile at 241 K. In addition to comparison to δ¹⁸O profiles for GRIP and other Greenland sites, the aim of the first simulation is to compare results from Crocus model to the models of Johnsen et al. (2000) and Bolzan and Pohjola (2000) run at this site with only diffusion along isotopic profiles. To compare our results to theirs, we consider an isothermal snowpack without meteorological forcing, and we deactivate modules of surface exchanges and heat transfer. The initial seasonal sinusoidal profile at GRIP is set using Eq. (19):

$$\begin{array}{}\text{(19)}& {\mathit{\delta }}^{\mathrm{18}}O(t,n)=-\mathrm{35.5}-\mathrm{8}\times \sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }\times z(t,n)}{a_{\text{GR}}\times {\mathit{\rho }}_{\text{ice}}/{\mathit{\rho }}_{\text{sn}}(t,n)}}\right),\end{array}$$

where z is the depth of the layer n, ρ_sn is its density, ρ_ice is the density of ice with a value of 917 kg m⁻³, and a_GR is the average accumulation at GRIP equal to 0.23 m ice eq. yr⁻¹ (Dahl-Jensen et al., 1993). The peak to peak amplitude value of 16 ‰ is close to the back-diffused amplitude at Summit (Sjolte et al., 2011).

The second simulation is run with evolving temperature in the snowpack. The snow temperature is computed by the model using meteorological forcing from ERA-Interim (see Table 4). In that case, the transport of isotopes in the vapor phase results from both diffusion along isotopic gradients and vapor concentration gradients. The initial snowpack is the same as in the previous simulation.

Figure 2 Simulation 1: attenuation of the seasonal δ¹⁸O_gcenter variation caused by diffusion along isotopic gradients in the vapor phase over 10 years (homogeneous and constant temperature of 241 K, original signal with a mean value of −35.5 ‰ and amplitude of 16 ‰). (a) Vertical homogeneous temperature profile; (b) δ¹⁸O profile at the beginning and end of the simulation; (c) deviation of the δ¹⁸O relative to the original profile for 10 dates; (d) evolution of the deviation to the original profile of δ¹⁸O.

In the two GRIP simulations, the modules of wind compaction and weight compaction are inactive. As weight compaction is taken to compensate for yearly accumulation (Eqs. 3 and 4), applying this compaction without precipitation would lead to an unrealistic drop in snow level. The wind compaction was absent from the model of Johnsen et al. (2000) and using this module would make comparisons more difficult.

3.3.2 Dome C simulations

In Simulations 3 to 6, we take advantage of the abundant documentation of the Dome C site to disentangle the different effects on the variations in water isotopic composition. All the simulations at Dome C were performed with an evolving temperature profile. Temperatures in the snow layers were computed using a modified meteorological forcing from ERA-Interim (Dee et al., 2011; Libois et al., 2014; see details in Table 4) and the modules of energy exchange and transfer. In this series of simulations, the δ¹⁸O values thus evolve as a result of diverging and/or alternating vapor fluxes. The simulations are ordered by increasing complexity. First, in Simulation 3, the modules of homogeneous compaction and wind drift are deactivated, as is the module of snowfall. Thus, the impact of vapor transport forced by temperature gradients on the snow isotopic compositions is clearly visible. Then, in Simulation 4, the module of compaction and the module of wind drift are activated to see their impact on the isotopes. We use an accumulation rate dm_sn∕dt for Dome C of 0.001 kg m⁻² per 15 min (see Eq. 3). Next, in Simulation 5, snowfall is added to assess how new layers affect snow δ¹⁸O values. Lastly, in Simulation 6, the model is run over 10 years at Dome C to build up a snowpack with realistic “sinusoidal” variation in δ¹⁸O values.

4 Results

4.1 Greenland

4.1.1 Results of the Crocus simulations (Simulations 1 and 2)

Figure 2 shows the result of Simulation 1, in which only diffusion along isotopic gradients is active, as in Johnsen et al. (2000). As expected the peak to peak amplitude of δ¹⁸O cycles is reduced because of diffusion. Over 10 years from 2000 to 2009, the amplitude decreases by 1.2 ‰, which corresponds to a 7.3 % variation.

Figure 3 shows the result of Simulation 2, i.e., with varying temperature in the snowpack. The attenuation is stronger than the one observed in the previous simulation. The minimum at 11.46 m increases by 1.03 ‰ over 10 years, and the maximum at 11.15 m decreases by 0.84 ‰. Thus, the total attenuation is ∼ 1.9 ‰ or 11.7 % for this height range. Below, the attenuation is smaller, with a total attenuation of only 6 % for heights between 10.54 and 10.85 m. If we compare attenuation for heights 11.46 and 11.56 m in the first and second simulation, we note that including temperature gradients leads to an increased attenuation by 50 %.

Figure 3 Simulation 2: attenuation of the seasonal δ¹⁸O_gcenter variation caused by diffusion in the vapor phase over 10 years (with temperature evolution, original signal with a mean value of −35.5 ‰ and amplitude of 16 ‰). (a) Vertical temperature profile for each summer; (b) δ¹⁸O_gcenter profile for each summer; (c) evolution of the deviation to the original profile of δ¹⁸O_gcenter. Note that temperature evolves during the whole year (see Fig. S1).

Figure 4 Evolution of the δ¹⁸O semi-amplitude with depth in shallow cores at NEEM, GRIP, and NGRIP (Steen-Larsen et al., 2011; Masson-Delmotte et al., 2015; White et al., 1997; Johnsen et al., 2000). The attenuation of the semi-amplitude values with depth was fitted using an exponential equation: $A=A_{\mathrm{0}}\exp \left(-\mathit{\gamma }z\right)-b$ with A₀=4.976 ‰, γ=0.08094, $b=-\mathrm{1.56}$ ‰ at GRIP and A₀=4.685 ‰, γ=0.06622, $b=-\mathrm{2.44}$ ‰ at NEEM. The dotted curve corresponds to the simulated attenuation at GRIP based on the Johnsen et al. model (diffusion length σ from their Fig. 2 and wavelength λ fitted on the Eurocore core from GRIP; White et al., 1997).

Between 11.46 and 11.56 m, the δ¹⁸O_gcenter values increase over 10 years by 1 to 4 ‰. This increase is not caused only by attenuation of the original sinusoidal signal. At h=11.60 m, the values get higher than the initial maximum, which was −36 ‰ at 11.64 m. There is therefore a local accumulation of heavy isotopes in this layer as a result of vapor transport. This maximum corresponds to a local maximum in temperature and is coherent with the departure of ¹⁸O-depleted water vapor from this layer. Thus, thermally induced vapor transport not only results in signal attenuation, but can also shift the δ¹⁸O value regardless of the initial sinusoidal variations.

Lastly, in the first 2–3 cm of the snowpack, strong depletion is observed over the period, with a decrease by 2 to 3 ‰ instead of 0.5 ‰ when the temperature gradients were absent (Simulation 1). This depletion probably results from the arrival of ¹⁸O-depleted water vapor from warmer layers below. This again shows the influence of temperature gradients that were absent from the previous simulation. However, note that in this simulation we neglect precipitation and the exchange of vapor with the atmosphere. Thus, the depletion observed here may not occur in natural settings when these processes are active.

In conclusion, at GRIP the diffusion of vapor as a result of temperature gradients has a double impact on isotopic compositions. It increases the attenuation in the first 60 cm of snow because of higher vapor fluxes. And it also creates local isotopic maxima and minima in a pattern corresponding to temperature gradients in the snowpack but disconnected from the original δ¹⁸O sinusoidal signal.

4.1.2 Comparison with core data

Here, we evaluate the attenuation of the initial seasonal signal in δ¹⁸O over 10 years at two Greenland ice-core sites, NEEM and GRIP. For the first site, we use four shallow cores (NEEM2010S2, NEEM2008S3, NEEM2007S3, NEEM2008S2) published in Steen-Larsen et al. (2011) and in Masson-Delmotte et al. (2015). For the second site, we use one shallow core (1989-S1) published in White et al. (1997). For the GRIP core, only the first 80 m is considered. Therefore, the data presented correspond to deposition and densification conditions like the modern ones. For NEEM the values of the four cores are taken together. For NEEM and GRIP, the semi-amplitude is computed along the core. In the first 10 m, the maximum value every 30 cm is retained, and deeper in the firn because of compaction, the maximum value every 20 cm is retained (see also the Supplement; Fig. 4). For this study, we have chosen to estimate attenuation in years with a clearly marked seasonal cycle, a strategy that can be debated but at least documented. Consequently, from this first series of maxima, a second series of maxima is computed with a larger window of 1 m. The “attenuated amplitudes” at each level are then defined as the ratio between these 1 m maxima and the initial 1 m maxima. Maximum semi-amplitudes every 5 m are also computed and displayed in Fig. 4. The 2.5 m attenuation is slightly higher at GRIP, leading to a remaining amplitude of 86 %, than at NEEM where the remaining amplitude is 90 % (Fig. 4). The amplitude decreases with depth in parallel for the two cores, with the amplitude at NEEM always staying higher than at GRIP. For comparison with our model, we estimate attenuation after 10 years, i.e., at a depth of ∼ 5.8 m for NEEM and ∼ 5.65 m for GRIP. The remaining amplitude is 80 and 72 % at GRIP and NEEM, respectively. Our Simulation 1 produced 7 % of attenuation only in the same duration, showing that our model run on an isothermal snowpack underestimates the attenuation observed in the data.

4.1.3 Comparison with other models

At 2.5 m at NGRIP, Johnsen et al. (2000) simulate a remaining amplitude of 77 % (Fig. 4). For a depth of 5.43 m corresponding to an age of 10 years, the simulated remaining amplitude is 57 %. For Bolzan and Pohjola (2000) at GRIP after 10 years, 70 % of the initial amplitude is still preserved. The slower attenuation for Bolzan and Pohjola (2000) compared to Johnsen et al. (2000) may be due more to the different sites considered than to the different models. GRIP has higher accumulation rates that should limit diffusion. Nevertheless, the attenuation of 30 % simulated by Bolzan and Pohjola (2000) at GRIP is stronger than the attenuation of 7 % simulated in our model. Town et al. (2008, Sect. 2.1) found attenuations of a few tenths per mil after several years when implementing only diffusion, a result consistent with ours since we get a decrease by 1.2 ‰ after 10 years.

We explore below the reasons for discrepancies between models. The equation for the effective diffusivity of vapor in firn used in our study is different from the ones used by Johnsen et al. (2000) and by Bolzan and Pohjola (2000). We do not consider the tortuosity factor l or the adjustable scale factor s of Bolzan and Pohjola (2000). However, using the values given by the previous authors for l and s leads to D_eff values ranging from $\mathrm{6.7}\times {\mathrm{10}}^{-\mathrm{6}}$ to $\mathrm{9.9}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ for a density of 350 kg m⁻³ and a temperature of 241 K, which is coherent with our value of $\mathrm{8.7}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹. As indicated by Bolzan and Pohjola (2000), the choice of one equation or another has little impact here.

The most probable difference lies in the way diffusion is taken into account. Johnsen et al. (2000) and Bolzan and Pohjola (2000) use a single equation of diffusion to predict the evolution of the isotopic composition of the layer. In our case, we specifically compute the fluxes in the vapor each second and at each depth level and deduce the evolution of δ¹⁸O in the grain center after sublimation–condensation and recrystallization. Denux (1996) and van der Wel et al. (2015) indicate that the model developed by Johnsen (1977) and used in Johnsen et al. (2000) overestimates the attenuation compared to observed values. For Denux (1996), the model of Johnsen (1977) should consider the presence of ice crusts and maybe also the temperature gradients in the surface snow to get closer to the real attenuation at remote Antarctic sites. Van der Wel et al. (2015) have compared the model results to a spike-layer experiment realized at Summit. Because an artificial snow layer cannot be representative of natural diffusion, they took care to evaluate diffusion based only on the natural layers present above and below the artificial layer. Van der Wel et al. (2015) propose three causes for the discrepancy between the Johnsen et al. model prediction and actual measured attenuation at GRIP. They blame either ice crusts, bad knowledge and parameterization of the tortuosity in the first meters of snow, and/or a bad description of the isotopic heterogeneity within the ice grain. In our model, the grain heterogeneity is included. Even if the parameters defining the mixing between the two compartments are not very well constrained (see Sect. 4.3), the attenuation is indeed smaller compared to the Johnsen model.

4.2 Dome C (Antarctica)

The aim of the Simulations 3 to 6, run at Dome C, is to isolate diffusion from other effects affecting water isotopic composition, i.e., wind drift and compaction.

Figure 5 Simulation 3: evolution of temperature and δ¹⁸O values from January to December 2001. (a) Temperature profiles for the first day of each month; (b) temperature evolution in the snowpack; (c) δ¹⁸O change in the grain surface compartment; (d) δ¹⁸O change in the grain center compartment. Here, “δ¹⁸O change” is defined as the difference between δ¹⁸O at t and at the beginning of the simulation for the selected layer.

Figure 6 Simulation 4: cumulative change in δ¹⁸O_gcenter values (vapor transport, compaction, and wind drift active).

Figure 7 Simulation 5: cumulative change in δ¹⁸O values at the grain center (relative to t₀) over 6 months. Simulation with snowfall with varying δ¹⁸O (function of T_air), vapor transport active, wind, and weight compaction active.

4.2.1 Simulation 3: without precipitation, without wind drift, and without homogeneous compaction

Figure 5 presents the results of temperature evolution (panels a and b) and δ¹⁸O evolution (panels c and d) for Simulation 3. The main changes in δ¹⁸O_gsurf and δ¹⁸O_gcenter occur in summer (Fig. 5c and d). On the one hand, the first 20 cm of snow tend to become ¹⁸O enriched by +0.2 ‰ for the grain center compartment. On the other hand, the first centimeter becomes depleted by 1.0 ‰ for grain center. This pattern is coherent with the temperature profiles for the summer period (Fig. 5a). Vapor moves out of the warmest layers and toward colder layers where it condensates. This causes an increase in δ¹⁸O in warm layers and a decrease in colder layers. This pattern is also confirmed by snow density changes (see Fig. S2).

During winter, the temperature generally decreases toward the surface (Fig. 5a). Vapor transport is thus reversed in the first 20 cm, but this only slightly reduces the dispersion of δ¹⁸O_gcenter values. On the first of August, the temperature at the surface temporarily increases to 235 K. This warm event strongly modifies the temperature profile in the snowpack and therefore the pattern of vapor transport. It is associated with an increase in δ¹⁸O values at the surface, which is particularly visible for the δ¹⁸O_gsurf values (Fig. 5c).

Thus, vapor transport can modify δ¹⁸O values in surface snow, even in the absence of precipitation or condensation from the atmosphere. This mechanism could explain the parallel evolution of surface snow isotopic composition and temperature described by Steen-Larsen et al. (2014) and Touzeau et al. (2016) between precipitation events.

4.2.2 Simulation 4: without precipitation, with wind drift, and with homogeneous compaction

Compaction and wind drift are not presumed to directly modify the δ¹⁸O values. However, the change in densities and layer thicknesses slightly modifies the temperature profile and the diffusivities. These processes thus could have an indirect impact on δ¹⁸O values. Figure 6 shows δ¹⁸O_gcenter changes that are reduced compared to the simulation without wind drift and compaction. This is coherent with a decrease in the density changes associated with vapor transport in the case with compaction (see Fig. S5).

4.2.3 Simulation 5: with precipitation, with wind drift, with homogeneous compaction

In Simulation 5, we add precipitation to wind and weight compaction effects. Both snowfall and wind compaction are responsible for irregular changes, respectively positive and negative, in the height of the snowpack (Fig. 7). In the new deposited layers, the δ¹⁸O_gcenter values reflect the δ¹⁸O values in the precipitation. They vary as expected from −40 ‰ on 31 December to −59 ‰ in July (Figs. 7, S8). The effect of vapor transport is visible only in “old” layers that were originally homogeneous in terms of δ¹⁸O. These old layers, which were reaching the surface in January, have been buried below the new layers and are found from 11 cm of depth downward in December.

4.2.4 Simulations 6: 10-year simulation at Dome C

Simulations 6 corresponds to a simulation run over 10 years at Dome C, with variable δ¹⁸O in the precipitation. Over these 10 years, about 1 m of snow is deposited. At the end of the simulation, the vertical profile of δ¹⁸O in the new layers has an average value of −49.7 ‰ and a semi-amplitude of 4.5 ‰ (Fig. 8). Here we take into account all the maxima and minima at a vertical resolution of 9 cm of fresh snow. Based on the atmospheric temperature variations only, the isotopic composition in the precipitation should vary around an average value of −53.2 ‰, with a semi-amplitude of 8.6 ‰. The main reason for this difference is the precipitation amounts: large precipitation events in winter are associated with relatively high δ¹⁸O values. The vertical resolution chosen for the model of 2.5 cm may also contribute to the decrease in the semi-amplitude. Light snowfall events do not result in the production of a new surface layer but are integrated into the old surface layer. As expected, the peak to peak amplitude of δ¹⁸O variations is then further reduced as a result of the two vapor diffusion processes and of associated vapor–solid exchanges. The effect of vapor transport is relatively small. To help with visualization, we selected four layers and displayed the evolution of δ¹⁸O in these layers over the years (Fig. 8d). The selected layers were deposited during winter 2000 and during summer seasons 2002, 2004, and 2006.

Figure 8 Simulation 6: evolution of δ¹⁸O_gcenter values as a result of snowfall and vapor transport over 10 years (compaction is inactive; merging between layers is allowed but limited). (a) Temperature profiles at mid-January for each year. (b) δ¹⁸O_gcenter profile at mid-January for each year. (c) Repartition of δ¹⁸O_gcenter values as a function of time and depth. (d) Evolution of δ¹⁸O_gcenter values after burial for four selected layers (deposited in winter 2000 and summer 2002, 2004, 2006). Note that we do not present the evolution of snow composition in the first year after deposition because the thin snow layers resulting from precipitation are becoming merged.

Figure 9 Test of the sensitivity of the model to the ratio of mass between grain surface compartments and total grain and to the interval of mixing between the two compartments (GRIP).

For the layer deposited during winter 2000, there is an increase in δ¹⁸O values of about +0.8 ‰ over 10 years. The slope is irregular, with the strongest increases occurring during summers between November and February when vapor transport is maximal. The slope is also stronger when the layer is still close to the surface, probably because of the stronger temperature gradients in the first centimeters of snow (Fig. 8a). For the layers deposited during the summers, the evolution of δ¹⁸O values is symmetric to that observed for winter 2000. Over 10 years, i.e., between 2000 and 2009, the δ¹⁸O amplitude thus decreases by about 1.6 ‰. This corresponds to a decrease of 18 % relative to the initial amplitude in the snow layers. This is higher than the 7 % attenuation modeled in Greenland for constant temperature and the 11.7 % attenuation observed when including diffusion caused by temperature gradients (Sect. 4.1). However, the comparison between the two sites is not straightforward, because of differences in temperature and accumulation counteracting each other. On the one hand, at GRIP the diffusion is forced by low vertical gradients of δ¹⁸O of the order of 0.24 ‰ cm⁻¹. These are much smaller than the typical δ¹⁸O gradients at Dome C, which are close to 1.10 ‰ cm⁻¹. On the other hand, the temperature of 241 K at GRIP is higher than the 220 K measured at Dome C, thus favoring diffusion.

4.3 Sensitivity tests for duration of recrystallization

We have shown above that the attenuation of the isotopic signal seems too small, at least for the GRIP site. In parallel, the parameters τ and Δt_gsurf/center of the model associated with grain renewal could only loosely be estimated, leading to uncertainty in the attenuation modeling. In this section, we perform some sensitivity tests to quantify how δ¹⁸O attenuation can be increased by exploring the uncertainty range in the renewal of the snow grain (Table 5). Indeed, the assumed values for the ratio between the grain surface and the total mass of the grain τ may have been underestimated or overestimated. The same is true for the periodicity of mixing between these two compartments Δt_surf/center.

The sensitivity tests are first designed for Greenland sites and run for 6 months, with initial amplitude of the sinusoidal δ¹⁸O signal of 16 ‰ and a fixed temperature of 241 K in all the layers (Fig. 9). First, we use a periodicity of mixing Δt_surf/center of 15 days and vary the value for the mass ratio τ: $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$, 5 $\times {\mathrm{10}}^{-\mathrm{4}}$, 3.3 $\times {\mathrm{10}}^{-\mathrm{2}}$. In practice, for Δt_surf/center=15 days, we realize mixing on day 2 and day 16 of each month. Second, we use the usual value of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ for τ and change the periodicity of the mixing to 2 days.

In the first case, where $\mathit{\tau }=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ and the mixing occurs every 15 days, the grain surface compartment is very small. Its original sinusoidal δ¹⁸O profile disappears in less than 1 day due to exchanges with vapor (not shown). The impact on the grain center is then very small with an increase in the first minimum by ∼ $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ ‰ over 6 months (Fig. 9a). In this case, the attenuation due to diffusion is even reduced compared to the results displayed above.

In the second case, where $\mathit{\tau }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ and the mixing occurs every 15 days, the grain surface compartment is larger and the attenuation is slower. Thus, in the grain surface compartment, half of the original amplitude still remains at the end of the simulation (not shown). The impact on the grain center compartment is clearly visible with an increase in the first minimum by $\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{2}}$ ‰ after 6 months (Fig. 9b).

In the third case, with $\mathit{\tau }=\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{2}}$, and mixing every 15 days, the attenuation of the sinusoidal signal in the grain surface compartment is only of 1 % because the grain surface compartment is very large. On opposite, attenuation in the grain center is quite large, i.e., the first minimum increases by $\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{2}}$ ‰ after 6 months (Fig. 9c).

In the fourth case, with $\mathit{\tau }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ and mixing every 2 days, the first minimum increases by $\mathrm{4.1}\times {\mathrm{10}}^{-\mathrm{2}}$ ‰ after 6 months for the grain center compartment (Fig. 9d). It is similar to the attenuation observed in the third case.

The results of these sensitivity tests suggest that the impact of vapor transfer on the grain center isotopic compositions is maximized when the grain surface compartment is large and/or refreshed often. They also clearly show that using a small grain surface compartment such as $\mathit{\tau }=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ drastically reduces the impact on the grain center isotopic values. However, our best estimates for τ and Δt_surf/center were not chosen randomly (see Sect. 3.1.3). Moreover, the use of $\mathit{\tau }=\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{2}}$ or Δt_surf/center=2 days leads to a near doubling of the δ¹⁸O attenuation (see above). This is not yet sufficient to explain the gap between our model output for isothermal simulation and the data. However, if this doubling is applicable to the case with temperature gradients, the attenuation obtained might reach the one observed in the data at GRIP.

At Dome C, sensitivity tests show that we can increase the attenuation by a factor of 3 by reducing the mixing time from 15 to 2 days (Fig. 10b–c). Similarly, if the ratio τ is put at $\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{2}}$ instead of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$, attenuation is more than doubled over 3 years (Fig. 10d–e). Thus, at Dome C the values of τ and Δt_surf/center seem to more strongly affect the attenuation obtained compared to GRIP. This greater sensitivity at Dome C could result from the influence of temperature gradients or from steeper δ¹⁸O gradients caused by the low accumulation. The average layer thickness of 2 cm in the first meter corresponds to ∼ 4 points per year at Dome C, but 35 points per year at GRIP.

Figure 10 Test of the sensitivity of the model to the ratio of mass between the surface and grain center compartments and to the interval of mixing between the two compartments (Dome C).

4.4 Additional missing processes

In the previous sections, we have seen that model outputs for GRIP generally lead to smaller attenuation than observed in ice cores. To improve the model compatibility with data, two kinds of approaches are possible. On the one hand, it would be useful to realize simulations adapted to on-site experiments, such as the one by van der Wel et al. (2015). This would allow for the verification of how diffusion can be improved in the model. For instance, previous studies have suggested that water vapor diffusivity within the snow porosity may be underestimated by a factor of 5 (Colbeck, 1983), but this is debated (Calonne et al., 2014). On the other hand, we also believe that other processes should probably be considered to explain the remaining attenuation. Ventilation is an additional process that has already been implemented in the snow water isotopic model of Town et al. (2008) and Neumann (2003). Because of strong porosity and sensitivity to surface wind and relief, ventilation is probably as important as diffusion in the top of the firn, even if diffusion is expected to be more effective at greater depths. For the Dome C simulation (Fig. 8), the slope d(δ¹⁸O)∕dz decreases slowly, indicating that diffusion remains almost as active at 60 cm than at 10 cm of depth. Neumann (2003) indicates that at Taylor Mouth the diffusion becomes the only process of vapor transport below 2 m of depth. For Dome C, for a temperature gradient of 3 ^∘C m⁻¹, we compute an average speed due to diffusion of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹. This is comparable to an air speed due to wind pumping of about $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ within the top meters of snow at WAIS (Buizert and Severinghaus, 2016). We conclude that, in as much as these results can be applied to Dome C, the two processes would have a comparable impact at this site in the first meters of snow. The next step for Crocus-iso development is thus to implement ventilation. Finally, we are also aware that in Antarctic central regions, the wind reworking of the snow has a strong effect in shaping the isotopic signal. A combination of stratigraphic noise and diffusion could indeed be responsible for creating isotopic cycles of non-climatic origin in the firn (Laepple et al., 2017). Wind reworking may also contribute to attenuation by mixing together several layers deposited during different seasons.

5 Conclusions and perspectives

Water vapor transport and water isotopes have been implemented in the Crocus snow model, enabling the depiction of the temporal δ¹⁸O variations in the top 50 cm of the snow in response to new precipitation, the evolution of the temperature gradient in the snow, and densification. The main process implemented here to explain post-deposition isotopic variations is diffusion. We have implemented two types of diffusion in the vapor phase: (1) water vapor diffusion along isotopic gradients and (2) thermally induced vapor diffusion. The vapor diffusion between layers was realized at the centimetric scale. The consequences of the two vapor diffusion processes on isotopes in the solid phase were investigated. The solid phase was modeled as snow grains divided into two sub-compartments: (1) a grain surface sub-compartment in equilibrium with interstitial water vapor and (2) an inner grain only exchanging slowly with the surface compartment. We parameterized the speed of diffusion through the renewal time of a snow grain and the proportion of the two snow grain compartments.

Our approach based on a detailed snow model makes it possible to investigate at a fine scale the various processes explaining the variations in density and δ¹⁸O in the firn. We look specifically at the effect of the evolution of the temperature gradient, new snow accumulation, and compaction events linked to wind drift. Over the first 30 cm, the snow density variations are mainly driven by compaction events linked to wind drift. Vapor transport and long-term compaction have secondary effects. Below 30 cm, wind-drift-driven compaction is no longer visible. Because of a strong temperature gradient and low density, water vapor transport will have a significant effect down to 60 cm. δ¹⁸O is primarily driven by variations in the δ¹⁸O of precipitation as expected. The seasonal variations are then attenuated by water vapor transport and diffusion along isotopic gradients, with an increase in these effects at higher temperatures, i.e., during summer periods.

From 10-year simulations of the Crocus-iso model both at GRIP in Greenland and Dome C in Antarctica, we have estimated the post-deposition attenuation of the annual δ¹⁸O signal in the snow to about 7–18 % through diffusion. This attenuation is smaller than the one obtained from isotopic data on shallow cores in Greenland, suggesting missing processes in the Crocus model when implementing water vapor. It is also significantly smaller than the diffusion implemented by Johnsen et al. (2000), but some studies have suggested that the Johnsen isotopic diffusivity is too strong (Denux, 1996; Van der Wel et al., 2015).

We see our study as a first step toward a complete post-deposition modeling of water isotope variations. Several other developments are foreseen in this model. First, wind pumping is currently not implemented in the Crocus model. This effect, implemented in the approach of Neumann (2003) and Town et al. (2008), is expected to have a contribution as large as the effect of diffusion for the post-deposition isotopic variations. Second, in low accumulation sites like Dome C, wind scouring probably has an important effect on the evolution of the δ¹⁸O signal at depth through a reworking of the top snow layers (Libois et al., 2014). This effect has not been considered here but could be implemented in the model in the next years. It could also play a role in the preservation of anomalously strong δ¹⁸O peaks at Dome C (Denux, 1996).

Other short-term developments concern the implementation of the exchange of water vapor with the atmosphere through hoar deposition. This is particularly timely since many recent studies have explored the parallel evolution of the isotopic composition of water vapor and surface snow during summer in both Greenland and Antarctica (Steen-Larsen et al., 2014; Ritter et al., 2016; Casado et al., 2016a, b). Similarly, the implementation of the ventilation of the snowpack in the model is important, since this effect is expected to significantly contribute to signal attenuation.

Another aspect is to look at the post-deposition d-excess and ¹⁷O-excess variations in snow pits. Recent studies have shown that the relationship between ¹⁷O-excess and δ¹⁸O is not the same when looking at precipitation samples and snow pit samples in East Antarctica (Touzeau et al., 2016). This observation questions the influence of diffusion within the snowpack on second-order parameters such as ¹⁷O-excess. Indeed, ¹⁷O-excess is strongly influenced by kinetic-diffusion-driven fractionation, which may be quantified by the implementation of ¹⁷O-excess in our Crocus-iso model.