Articles | Volume 13, issue 9
Model description paper
17 Sep 2020
Model description paper |  | 17 Sep 2020

The Community Firn Model (CFM) v1.0

C. Max Stevens, Vincent Verjans, Jessica M. D. Lundin, Emma C. Kahle, Annika N. Horlings, Brita I. Horlings, and Edwin D. Waddington

Models that simulate the evolution of polar firn are important for several applications in glaciology, including converting ice-sheet elevation change measurements to mass change and interpreting climate records in ice cores. We have developed the Community Firn Model (CFM), an open-source, modular model framework designed to simulate numerous physical processes in firn. The modules include firn densification, heat transport, meltwater percolation and refreezing, water isotope diffusion, and firn-air diffusion. The CFM is designed so that new modules can be added with ease. In this paper, we first describe the CFM and its modules. We then demonstrate the CFM's usefulness in two model applications that utilize two of its novel aspects. The CFM currently has the ability to run any of 13 previously published firn densification models, and in the first application we compare those models' results when they are forced with regional climate model outputs for Summit, Greenland. The results show that the models do not agree well (spread greater than 10 %) when predicting depth-integrated porosity, firn age, or the trend in surface elevation change. In the second application, we show that the CFM's coupled firn-air and firn densification models can simulate noble gas records from an ice core better than a firn-air model alone.

1 Introduction

Snow that falls on an ice sheet transitions to ice through an intermediate stage called firn. Knowledge of the physics of firn densification has several applications in glaciology. Studies of ice-sheet mass balance using altimetry methods require knowing the mass, and mass changes, of the firn to estimate the contribution of the ice sheets to sea-level rise (Shepherd et al.2012; The IMBIE Team2018). Ice-core studies require knowledge of the age of the firn at the depth at which bubbles of air become trapped in order to determine the difference between the age of air in the bubbles and the ice that encloses the bubbles (called Δage; Blunier and Schwander2000). Both of these applications require a firn densification model. In addition, ice-core researchers use firn-air models to simulate the diffusion of atmospheric gases through the porous firn; these models can be used, for example, to estimate the age of gases when they become trapped in bubbles (e.g., Buizert et al.2012).

Firn is commonly divided by density into three zones based on the dominant physics of densification (Herron and Langway1980; Maeno and Ebinuma1983). The first zone extends from the surface, where density is often assumed to be  300–350, to 550 kg m−3. In zone 1, densification is usually considered to be due to grain-boundary sliding and settling (Alley1987). In zone 2, which spans the densities between 550 and  830 kg m−3, densification occurs due to sintering processes (Gow1975). Near a density of 830kg m−3, bubbles of air become trapped in the ice matrix. This density is referred to as the bubble close-off (BCO) density, which is reached at the corresponding BCO depth. The BCO depth varies by site. Further densification occurs due to compression of the bubbles in zone 3, which comprises the bubbly ice between the BCO density and the ice density. The densification rate slows significantly in zone 3 as the pressure in the bubbles increases (Goujon et al.2003).

Numerous models have been developed to describe the physics of firn densification. In addition to predicting the evolution of firn density, most firn densification models also simulate the firn's temperature evolution by coupling a heat diffusion model. A common way to form a firn densification model is to assume that, for a given site, the accumulation rate is constant and the firn-density profile is in steady state. Using this steady-state assumption, known as Sorge's law, the change in density ρ with depth (dρ∕dz) can be converted to a material-following (Lagrangian) change in density with time (dρ∕dt) (Bader1954). Using depth–density data from many sites, firn densification models can be formulated as a function of temperature (often through an Arrhenius term with a tuned activation energy), accumulation rate (a proxy for stress), and one or more tuning parameters.

Several firn densification models have been developed without invoking Sorge's law. Some of these have used firn strain-rate data; these include work by Arthern et al. (2010), who measured firn compaction rates in real time in Antarctica using “coffee-can” type strain gauges (Hulbe and Whillans1994), and by Morris and Wingham (2014), who inferred firn compaction rates by tracking layering in repeated high-resolution density logs of boreholes. Other studies have worked to develop a model based on the microphysical processes driving densification; e.g., Alley (1987) developed a model for zone 1 densification by applying grain-boundary sliding theory. Arnaud et al. (2000) combined the work by Alley (1987) with theory describing pressure sintering of spherical powders (Arzt1982) to simulate densification in zone 2.

The evolution of firn density is governed by grain-scale (i.e., microstructural) processes (Arnaud et al.2000; Morris and Wingham2014), but most firn models predict density evolution based only upon the accumulation rate and temperature. These fields can come from a regional climate model (RCM), e.g., the Regional Atmospheric Climate Model (RACMO;  Noël et al.2018) or Modèle Atmosphérique Régional (MAR;  Fettweis et al.2017), in situ weather stations such as the Greenland climate network (GC-Net;  Steffen and Box2001; Vandecrux et al.2018), or ice-core data (e.g.,  Buizert et al.2015). Firn densification models also need a surface-density boundary condition; this can be assumed to be constant through time (Fausto et al.2018), or it can be predicted, for example, by using an empirical parameterization based on temperature (e.g.,  Kuipers Munneke et al.2015) or some other variable.

Atmospheric gases move through the firn's pore space above the BCO depth. Below the BCO depth these gases are trapped in bubbles, and they preserve a record of past atmospheric composition. Gas transport in firn is commonly modeled by dividing the firn into three zones by dominant transport mechanisms, which differ from the three zones of densification (Sowers et al.1992). Near the surface is the convective zone, which may be from 0 to  20 m thick; in this zone, convective mixing due to wind pumping and buoyancy dominates gas transport and keeps the air well-mixed and in equilibrium with the free atmosphere above (Kawamura et al.2006). Below the convective zone is the diffusive zone, where gas transport is driven primarily by diffusion along chemical concentration gradients. Additionally, isotopic fractionation occurs in the diffusive zone due to gravitational and thermal effects: gravitational fractionation causes heavier isotopes to become enriched (i.e., their relative abundance increases) at greater depths, and thermal fractionation causes heavier isotopes to become enriched at the cold end of a temperature gradient (Severinghaus et al.1998). Numerous models have been developed to simulate air movement in firn; they aid in using firn-air measurements to reconstruct past atmospheric conditions (Buizert et al.2012).

We have developed the Community Firn Model (CFM), an open-source model framework that includes a suite of published firn densification models, a firn-air model, and numerous modules to simulate other physical processes in firn. We created the CFM to be a resource to the glaciological community at large. We recognize that many research groups have their own firn models, but our goals with the CFM are (1) to provide a model to research groups who need firn model outputs but do not want to code a model themselves, (2) to provide a point of reference for research groups to compare their model output against, and (3) to enable firn densification model comparisons within a single model framework to improve understanding of firn model uncertainties within various applications. In this paper, we describe the model and demonstrate its utility in two model applications.

2 The Community Firn Model

The CFM is an open-source, modular firn model framework. It is coded in Python 3 and is available for download on GitHub. It is designed to simulate numerous processes associated with firn; “modular” refers to the fact that the CFM was constructed so that users can choose which of these processes they would like to simulate in a given model run, and a module is a piece of code that simulates a particular process (e.g., density evolution). Modules to simulate additional processes can be added with minimal alteration of existing code. The CFM's modularity allows the user to easily choose which physical processes to simulate and which model outputs to save in a particular CFM run. The core modules of the CFM track the evolution of the firn density and temperature. Other modules simulate grain-size evolution, firn-air diffusion, water isotope diffusion, meltwater percolation and refreezing, and layer thinning due to horizontal strain; these latter modules also require execution of the core modules.

2.1 CFM workflow

Prior to running the CFM, the user sets the parameters specific to the model run in a JSON-formatted configuration file. These include, among others, which firn densification physics to use, the time-step size, and the model domain thickness. The configuration file also includes the paths to the files used for forcing the model (i.e., the surface boundary conditions). The CFM's GitHub repository includes an example configuration file preset with default values, and the CFM's documentation includes detailed descriptions of each of the parameters.

The CFM is forced at the upper boundary (i.e., the ice-sheet surface) using surface temperature (i.e., the temperature of the snow at the surface), surface density, accumulation-rate data, and any other surface boundary condition needed for a particular module. These fields are input via a CSV file that includes time in the first row and value in the second row. If the times in the input files are not the same (e.g., as might be the case with climate data from ice cores), the CFM interpolates them onto a common timeline. Many RCM outputs, including RACMO and MAR, are stored in NetCDF files. We include a script on the CFM's GitHub repository to assist users in creating CSV climate-forcing files from NetCDF data, and a goal for future releases is to allow users to run the CFM with direct input from the NetCDF files.

The CFM uses a Lagrangian grid with a fixed number of model volumes; each volume represents a layer of firn with uniform properties. The number of volumes is determined by the thickness of the model domain, the time-step size, and the mean annual accumulation rate; one volume is the accumulation for a single time step. At each time step, accumulation is added as a new volume at the surface, and a volume is removed from the bottom of the grid. There are no limitations on the time-step size or thickness of the model domain. However, because the thickness of a model volume depends directly on the amount of accumulation during a time step, the number of volumes increases as the time-step size decreases and increases the computational burden. The CFM has an optional scheme to merge model volumes and thus reduce computing time.

A CFM run begins by first “spinning up” the model to a steady state, which becomes the initial condition for the “main” model run. A CFM run begins by setting the depth–density profile to that predicted by the Herron and Langway (1980) analytic firn densification equations (Sect. 2.2.1) using the site's mean accumulation and temperature. The spin-up then evolves the firn by time stepping forward using the specified firn densification equation. The climate forcing during spin-up can either be constant (appropriate for runs with large time steps, e.g., runs to simulate Δage through time for ice cores) or can include climatological noise (appropriate for runs with small time steps, e.g., simulating surface elevation change). Although the user specifies how long the spin-up should last, it is recommended to spin up long enough to reset the entire grid (i.e., to flush out all of the initial volumes and replace them with new volumes). During spin-up, the CFM evolves the density, temperature, age, and other properties that might be included in a model run (e.g., grain size).

After the spin-up has completed, the main model run begins. It operates in the exact same way as the spin-up, except the model is forced with varying temperature, accumulation rate, and other boundary conditions, rather than the steady-state values. When the model run is complete, the model outputs are saved in a single HDF5-formatted file. The user specifies which model outputs to save; the options are firn depth, density, age, temperature, compaction rate, grain size, water isotope values, BCO depth and age, depth-integrated porosity, liquid water content, and gas concentrations. The resolution of the model outputs is specified by the user in the JSON configuration file; by default the CFM saves the outputs on the entire model grid at each time step.

We next describe the various modules built into the CFM, with particular focus on the firn-density and firn-air modules.

2.2 Density

The CFM is coded to include 13 previously published firn densification equations (listed in Table 1), and it is designed so that it is easy for the user to choose which firn densification equation to use in a particular CFM run. We note that the word “model” is often used to refer to an equation, which can lead to ambiguity. We use “CFM” to refer to the entire Community Firn Model framework, and we use “firn densification model” to refer to an equation or set of equations that simulates the physics of firn densification. Thus, running the CFM includes implementing a firn densification model.

Herron and Langway (1980)Barnola et al. (1991)Goujon et al. (2003)Li and Zwally (2011, 2015)Helsen et al. (2008)Arthern et al. (2010)Ligtenberg et al. (2011)Kuipers Munneke et al. (2015)Simonsen et al. (2013)(Vionnet et al.2012)Morris and Wingham (2014)

Table 1List of firn densification models and their abbreviations coded in the CFM and included in the study detailed in Sect. 3. Arthern et al. (2010) describe two different models; see Sect. 2.2.6 for details.

Download Print Version | Download XLSX

A general form used in many firn densification models assumes that a firn layer's change in density ρ through time t is a function of the temperature T, accumulation rate b˙, and current density:

(1) d ρ d t = f ( T , b ˙ , ρ ) .

Density evolution in the CFM is handled with an explicit numeric scheme, i.e.,

(2) ρ new = ρ old + ( d ρ / d t ) d t .

Most of the firn densification models in the CFM use the accumulation rate as a proxy for the stress. If the accumulation rate is constant in time, the overburden stress σ at depth z is related to the mass accumulation rate b˙ by the relation σ(z)=b˙gτ(z), where g is gravity and τ(z) is the age of the firn at depth z. For the models that are forced with the accumulation rate (as opposed to stress), the CFM by default uses the mean accumulation rate b˙ over the lifetime of each parcel of firn rather than the instantaneous accumulation rate at a given time step. The mean accumulation rate for any layer of firn at time t and depth z with agez,t is determined by the integrated accumulation-rate (b˙) history (Li and Zwally2011, 2015):

(3) b ˙ ( z , t ) = 1 age z , t t - age z , t t b ˙ ( t ) d t .

The CFM uses b˙ because a firn densification model dependent on the instantaneous accumulation rate will predict that no densification occurs when the accumulation rate is zero (Li and Zwally2011), which is not realistic. This approach may be different than how some of the models were originally formulated, and the CFM includes an option to use the instantaneous accumulation rate. In steady state, the mean accumulation rate is the same as the instantaneous rate.

The surface density ρs of a new layer of firn in the CFM can be a constant value or can vary in time. In the case of time-varying ρs, it is determined by a parameterization (e.g.,  Kuipers Munneke et al.2015) or by randomly selecting a value from a specified distribution.

We have coded each of the firn densification models in the CFM as we have interpreted their descriptions in their original publications, and we have corrected any known errors. We next provide a basic description of each of the models and any nuances associated with coding them. The subsection headings also include the abbreviations that we use for each model in the application described in Sect. 3. For additional descriptions of the firn densification models included in the CFM, see the original publications.

2.2.1 Herron and Langway (1980, HL)

Herron and Langway (1980) is a benchmark firn densification model (Lundin et al.2017); nearly all firn densification models developed since 1980 are based in part on assumptions made by those authors. They used Sorge's law and depth–density data from 17 firn cores to derive a widely applicable firn densification rate equation. The CFM includes three formulations of the Herron and Langway (1980) model, which are detailed in Lundin et al. (2017): a “dynamic” model, a “stress-based model”, and an “analytic” model. In steady state, all three give the same result, but the outputs vary in transient simulations. The CFM uses the analytic model to generate an initial condition. Here we describe only the “dynamic model”, which is used in the application in Sect. 3.

Two assumptions in the Herron and Langway (1980) model have been used in numerous other firn densification models. They are the following: (1) the change in porosity is linearly related to the stress change resulting from new snow accumulation (i.e., the densification rate is a function of the porosity; Schytt1958; Robin1958), and (2) the firn's densification rate has an Arrhenius dependence on temperature. These assumptions can be incorporated into a densification-rate equation:

(4) d ρ d t = c ( ρ ice - ρ ) ,


(5) c = k exp - Q R T b ˙ a ,

where k and a are constants, Q is the Arrhenius activation energy (kJ mol−1), R is the gas constant (8.314 kJmol-1K-1), and T is the temperature (K). For HL, c in Eq. (4) is given by

(6) c = c 0 = 11 exp - 10.16 R T b ˙ 1.0 ( ρ 550 kg m - 3 ) = c 1 = 575 exp - 21.4 R T b ˙ 0.5 ( ρ > 550 kg m - 3 ) .

HL uses units of meters water equivalent per year (m w.e. a−1) for b˙. T in the original model was the mean annual site temperature Tm, but since its establishment HL has also been implemented such that T is the temperature of a specific parcel of firn. The parameters k, a, and Q were all tuned to best fit the firn-core data. The activation energy derived by Herron and Langway (1980) is lower than most other models, which causes it to be less sensitive to sub-annual temperature variability. We note that because of the different values of the exponent a on b˙, k in Eq. (5) has different units for zones 1 and 2. The units for ρ in HL are megagrams per cubic meter (Mg m−3), which are numerically equivalent to units grams per cubic centimeter (g cm−3).

2.2.2 Barnola et al. (1991, BAR)

The Barnola et al. (1991) model was developed for ice-core Δage calculations. It uses the Herron and Langway (1980) model for zone 1 densification. For ρ>550kg m−3, the densification rate is given by

(7) d ρ d t = ρ i A 0 exp - Q R T f σ d e f f n ,

where A0=2.54×104MPa-3s-1, the activation energy Q is 60kJ mol−1, σeff is the effective stress (MPa), and n=3. For zone 2 densification, Barnola et al. (1991) derived an equation empirically to match the densification rate and its derivative at the zone 1 to zone 2 and bubble close-off transitions, and f is given by

(8) f = 10 α ρ 3 + β ρ 2 + δ ρ + γ ( 550 ρ 800 kg m - 3 ) ,

with α=-37.455, β=99.743, δ=-95.027, and γ=30.673.

In zone 3, beyond the 800 kg m−3 density horizon, f is taken from Pimienta (1987):

(9) f = 3 16 ( 1 - ρ / ρ i ) / ( 1 - ( 1 - ρ / ρ i ) 1 / 3 ) 3 ( ρ > 800 kg m - 3 ) .

2.2.3 Arnaud et al. (2000) and Goujon et al. (2003, GOU)

Arnaud et al. (2000) developed a densification model based upon descriptions of grain-scale physical processes in firn, and Goujon et al. (2003) extended that model by adding a heat diffusion component. These models were developed for ice-core Δage reconstructions in Antarctica. In the CFM we refer to this family of models as the Goujon model because the CFM includes a heat diffusion module. For zone 1 densification, the Goujon et al. (2003) model is based on grain-boundary sliding work by Alley (1987). It describes the densification rate in zone 1 as

(10) d D d t = γ P D 2 1 - 5 3 D ( ρ 550 kg m - 3 ) ,

where D=ρ/ρice is the relative density, P is the overburden pressure (bar), and γ is a scaling factor that depends on the viscosity of grain boundaries and the geometry of the grains. It is notable that GOU's densification rate in zone 1, unlike most other firn densification models, does not depend on temperature.

Goujon et al. (2003) base their description of zone 2 densification on sintering theory from Arzt (1982):

(11) d D d t = 4.1817 × 10 4 exp - E A R T D 2 D 0 1 / 3 a π 1 / 2 4 π P 3 a Z D 3 ( ρ > 550 kg m - 3 ) ,

where EA is the activation energy given as 60 kJ mol−1, R is the gas constant, T is the temperature (K), a is the average contact area between the grains relative to the initial grain radius, and Z is the coordination number, i.e., the average number of neighboring grains to a central grain in the firn crystalline structure. D0 is the zone 1–zone 2 transition relative density. Details on determining a and Z can be found in the original publications. Unlike other firn densification models, which specify a constant transition density, the Goujon et al. (2003) model uses a transition density that depends on Tm (K), given by

(12) D 0 = 0.00226 T m + 0.03 .

Goujon et al. (2003) specify that γ in Eq. (10) should be set so that the densification rate is continuous at D0.

Buizert et al. (2015) described an issue in implementing the Goujon et al. (2003) model, which we review here. In the event that D0≥0.6, dD∕dt given by Eq. (10) becomes zero for D=0.6 and negative for D≥0.6, which is not realistic. Additionally, at D=D0, the densification rate predicted by Eq. (11) is infinite because the contact area a equals zero. We avoid these issues in our implementation of GOU by doing the following. We limit the value of D0 given by Eq. (12) to a maximum value of 0.59, which occurs for temperatures greater than -25C. This value of D0 corresponds to a density of 541kg m−3, which results in GOU always predicting the zone 1–zone 2 transition occurring at lower densities than the commonly used value of 550 kg m−3. We follow the suggestion in Buizert et al. (2015) and put the zone 1–zone 2 transition at relative density D0=D0+ϵ, where ϵ is a small number. The densification rate in zone 2 is still calculated using Eq. (11) using D0 given by Eq. (12). We then iterate to find γ in Eq. (10) that gives the maximum dD∕dt at the bottom of zone 1, with the condition that it does not exceed dD∕dt given by Eq. (11) at the top of zone 2.

2.2.4 Li and Zwally (2011, LZ11) and Li and Zwally (2015, LZ15)

The Li and Zwally (2011) and Li and Zwally (2015) models are the latest in a lineage of models developed by the authors (Li and Zwally2002, 2004). The models were developed to predict the surface elevation changes associated with seasonal variability in accumulation and firn compaction rates. LZ11 and LZ15 are tuned to model firn in Greenland and Antarctica, respectively. Both share the same basic form:

(13) d ρ d t = β 8.36 ( 273.2 - T K ) - 2.061 b ˙ ( ρ i - ρ ) ,

where TK is the firn temperature as a function of time and depth in Kelvin, and b˙ is the mean accumulation rate over the lifetime of a parcel of firn (Eq. 3) (m w.e. a−1). The difference between LZ11 and LZ15 is in the parameter β. For LZ11,

(14) β = β 1 = - 9.788 + 8.996 b ˙ m - 0.6165 T m,c ( ρ 550 kg m - 3 ) = β 2 = β 1 / ( - 2.0178 + 8.4043 b ˙ m - 0.0932 T m,c ) ( ρ > 550 kg m - 3 ) .

For LZ15,

(15) β = β 1 = - 1.218 - 0.403 T m,c ( ρ 550 kg m - 3 ) = β 2 = β 1 ( 0.792 - 1.080 b ˙ m + 0.00465 T m,c ) ( ρ > 550 kg m - 3 ) ,

where Tm,c is the mean annual surface temperature in Celsius and b˙m is the long-term accumulation rate at the site being modeled.

LZ11 and LZ15 predict unrealistically high densification rates for firn near the freezing temperature (and infinite at the freezing temperature), which makes them unsuitable for simulations of wet firn.

2.2.5 Helsen et al. (2008, HEL)

The Helsen et al. (2008) model was developed to simulate firn-column thickness changes in Antarctica to improve ice-sheet mass change estimates derived from satellite altimetry observations. Its development was based on the work of Li and Zwally (2002) and uses the same general form for its densification equation (Eq. 13). Helsen et al. (2008) used additional firn-core data from Antarctica to derive a different value for β. HEL uses a single β for zone 1 and zone 2 densification, given by

(16) β = β 1 = β 2 = 76.138 - 0.28965 T m .

The mean annual surface temperature Tm in Eq. (16) is in Kelvin.

2.2.6 Arthern et al. (2010, ART-T & ART-S)

Arthern et al. (2010) derived a firn densification model using firn compaction rate data from several sites in the Filchner–Ronne sector of Antarctica. The model is notable because it was the first firn densification model to be based on compaction-rate measurements rather than upon a derived compaction rate from Sorge's law. The authors also identified two processes in firn, diffusion of water molecules through the ice lattice and grain growth, that have different activation energies. They hypothesized that these processes acting in concert result in a lower effective activation energy, which could explain the low temperature sensitivity in HL.

Arthern et al. (2010) describe two implementations of their model: the first is the complete transient, dynamical model described in their appendix, which we refer to as ART-T. It includes equations for densification rate based on evolving stress σ and grain radius r (based on the work of Gow et al.2004):

(17) d ρ d t = k c ( ρ i - ρ ) exp ( - E c / R T ) σ / r 2 d r 2 d t = k g exp ( - E g / R T ) d σ d t = b ˙ g ,

where kc and kg are empirically derived constants. The water-molecule diffusion activation energy Ec is 60 kJ mol−1, and the grain-growth activation energy Eg is 42.4 kJ mol−1. ART-T is not in common use; it is sensitive to the surface grain size, which is poorly constrained for model simulations over the range of climates encountered on ice-sheet scales.

Arthern et al. (2010) use several simplifying assumptions, including that of steady accumulation, to derive the second implementation of their model, which is the model presented in their main text. We refer to this implementation as ART-S. The densification equations for zone 1 and zone 2 use the same form as the Herron and Langway densification equation (Eq. 4), with parameters c given by

(18) c 0 = 0.07 b ˙ g exp - E c R T + E g R T m ( ρ 550 kg m - 3 ) c 1 = 0.03 b ˙ g exp - E c R T + E g R T m , ( ρ > 550 kg m - 3 ) ,

with b˙ the mass accumulation rate (kgm-2a-1) and Tm in Kelvin. ART-S forms the basis of the models described by Ligtenberg et al. (2011), Simonsen et al. (2013), and Kuipers Munneke et al. (2015).

2.2.7 Ligtenberg et al. (2011, LIG) and Kuipers Munneke et al. (2015, KM)

The Ligtenberg et al. (2011) and Kuipers Munneke et al. (2015) models comprise the subsurface scheme of the regional climate model RACMO (Noël et al.2018; van Wessem et al.2018). They were developed to simulate firn densification and meltwater percolation and refreezing in Antarctica and Greenland, respectively. The development of their densification equations was based on ART-S, but the authors used firn-core data to widen the ART-S applicability across the ice sheets. LIG and KM use the same form as ART-S with the exception that c0 and c1 in Eq. (18) are multiplied by additional tuning coefficients. For LIG,

(19) c 0 (LIG) = [ 1.435 - 0.151 ln ( b ˙ ) ] c 0 (ART-S) c 1 (LIG) = [ 2.366 - 0.293 ln ( b ˙ ) ] c 1 (ART-S) ,

and for KM,

(20) c 0 (KM) = [ 1.042 - 0.0916 ln ( b ˙ ) ] c 0 (ART-S) c 1 (KM) = [ 1.734 - 0.2039 ln ( b ˙ ) ] c 1 (ART-S) .

For LIG and KM, b˙ has units of kilograms per square meter per year (kgm-2a-1). Ligtenberg et al. (2011) specified a densification rate as a function of both the firn temperature and mean annual surface temperature, as is done in ART-S (Eq. 18). Kuipers Munneke et al. (2015) used T rather than Tm in the grain-growth term of the Arrhenius factor; Steger et al. (2017) modified this to use Tm, and we include this latest version in the CFM.

2.2.8 Simonsen et al. (2013, SIM)

The Simonsen et al. (2013) model was developed for ice-sheet mass balance studies. The authors used a Monte Carlo inverse method with airborne radar data and regional climate model data to tune the parameters in the firn densification model. Like LIG and KM, SIM also uses the densification equation from ART-S as a basis and multiplies c0 and c1 in Eq. (18) by additional factors:

(21) c 0 (SIM) = f 0 c 0 (ART-S) c 1 (SIM) = f 1 61.7 b ˙ 0.5 exp - 3800 R T m c 1 (ART-S) .

Simonsen et al. (2013) found that the parameters f0 and f1 needed to be tuned and/or specified to simulate the firn at a particular site (see  Lundin et al.2017, Appendix); the CFM uses default values f0=0.8 and f1=1.25 (Sebastian Simonsen, personal communication, 2015).

2.2.9 CROCUS (Brun et al., 1992; Vionnet et al., 2012), CRO

The Crocus model was developed for mountain snowpacks (Brun et al.1992), but it has also been used to simulate firn densification and hydrology (e.g.,  Langen et al.2017; Verjans et al.2019). Its equations are also used for the subsurface scheme in the RCM MAR (Cullather et al.2016), which is used to simulate the surface mass balance of the Greenland and Antarctic ice sheets (Fettweis et al.2017; Agosta et al.2019; Alexander et al.2019).

The Crocus model gives compaction in terms of a constitutive equation, relating stress σ to the densification rate with a viscosity η (Vionnet et al.2012):

(22) d ρ d t = ρ σ η η = f 1 f 2 η 0 ρ c η exp a η ( T melt - T ) + b η ρ ,

where ρ is the density, Tmelt is the melting temperature, η0=7.62237×106kgs-1m-1, aη=0.1K-1, bη=0.023m3kg-1. f1 and f2 are snow-viscosity correction factors. f1 accounts for viscosity differences due to the presence of liquid water; it is set to 1 in the CFM for model simulations in the dry-firn zone. f2 accounts for angular grains. Following Langen et al. (2017) and van Kampenhout et al. (2017), f2 is set to 4 in the CFM. Vionnet and others (2012) give cη=250kgm-3, but also following van Kampenhout et al. (2017), we set cη=358kgm-3 in the CFM.

2.2.10 Morris and Wingham (2014, MW)

Morris and Wingham (2014) used a neutron probe to measure firn density in boreholes in successive years at several sites on the Greenland Ice Sheet. The high vertical spatial resolution of these measurements allowed the authors to measure both the strain between firn layers and density changes in those layers. They used this information to derive a compaction equation for dry firn of density less than 550 kg m−3:

(23) ϵ ˙ = k 0 * ρ w g ρ i - ρ ρ 1 - M 0 m 1 H ( τ ) exp - E α R T σ ,

where H(τ) is a “temperature history function”:

(24) H ( τ ) = τ 0 τ exp - E α R T ( τ ) d τ .

In the above equations, ρi and ρw are the densities of ice and water, respectively, k0* is a densification constant, and τ is the age of a parcel of firn. τ0 is the age of the firn after it leaves the thin surface snow layer, which we take in the CFM to be zero. m is the normalized deviation of ρ(z) from a quadratic curve fit to the density profile, and M0 is a scaling constant. The CFM by default uses the preferred activation energy Eα presented in Morris and Wingham (2014) of 110 kg m−3, though the authors also test their model with values of 60 and 200 kJ mol−1. As such, the CFM is coded to allow the user to vary Eα easily.

Morris and Wingham (2014) described simplifying assumptions for their model, which allowed them to use the same densification constant k0*=11mw.e.-1 as Herron and Langway (1980). However, this description was based on an error in the calculation of the history function from their data. Conceptually, the result of the error is that the authors' original hypothesis of a single densification process does not hold. In practice, this error causes the model to predict unrealistic densification rates. Table 2 shows the corrected values of EαEH (revising Table 2 from Morris and Wingham2014), which can be used to calculate the corrected value of k0* with the following equation:

(25) k 0 * = 11 exp - E * - ( E α - E H ) R T m ,

where E* is the Herron and Langway (1980) activation energy.

Table 2EαEH (kJ mol−1), revised from Table 2 of Morris and Wingham (2014).

Download Print Version | Download XLSX

The error and its correction were described in a personal communication with Elizabeth Morris in April 2019; that communication is included with the CFM's documentation. MW only specifies densification rates for zone 1; the CFM is coded to use the Herron and Langway (1980) equation for zone 2 densification.

2.3 Temperature evolution

The temperature in firn evolves by diffusion and advection. Heat diffusion in the CFM is modeled using a fully implicit finite-volume scheme (Patankar1980). Advection of heat is inherently handled by the Lagrangian scheme. The CFM uses a Dirichlet (prescribed temperature) boundary condition at the surface; the surface temperature is set by the input forcing data. At the bottom of the domain the CFM uses a Neumann (i.e., prescribed gradient, set to zero) boundary condition for the temperature by default; this condition can be easily adapted to use a nonzero gradient or a Dirichlet boundary condition.

Simulating heat diffusion in firn requires knowledge of the thermal conductivity, but there is no universally agreed upon parameterization for the thermal conductivity of snow and firn. The CFM includes a number of parameterizations for the thermal conductivity that have been published previously. Those are the following: Anderson (1976), Brandt and Warren (1997), Jiawen et al. (1991), Lüthi and Funk (2001), Riche and Schneebeli (2013), Schwander et al. (1997), Schwerdtfeger (1963), Sturm et al. (1997), Van Dusen (1929), and Yen (1981).

In the wet-firn zone of an ice sheet, there is additional heat transport due to advection of liquid water at 0 C through the firn's pore space and due to latent heat release from meltwater refreezing. The CFM simulates the advective component with a meltwater percolation scheme, which is described in Sect. 2.6. The latent heat from refreezing is handled in one of two ways, depending on the meltwater percolation scheme that is used. The first uses a fully implicit, finite-volume, enthalpy diffusion scheme to resolve latent heat release and heat diffusion (Voller et al.1990). We note that a similar enthalpy-based method was employed by Meyer and Hewitt (2017). The second computes latent heat release in the meltwater percolation scheme and separately uses the dry-firn heat diffusion scheme; details are provided in Verjans et al. (2019).

The CFM does not incorporate a scheme to account for the impact of shortwave radiation penetration into the firn, although research has suggested that it can affect the temperature in the near-surface snow by several degrees (Kuipers Munneke et al.2009). Adding a module to account for this effect could be an area for future development of the CFM. This would require the implementation of a scheme that computes the transfer of these radiative components into firn, forced by surface values that must be provided by regional climate models or weather station data.

2.4 Water isotope diffusion

The CFM includes a module that calculates the diffusion of water isotopes, which occurs due to sublimation and deposition of water molecules in the firn. This process is important to the interpretation of ice-core records. Variability in the water isotopic composition of snow crystals that fall on the surface is damped by a diffusion process as the snow advects downward through the firn column. In the vapor phase, water molecules diffuse through the firn column's interconnected pore space, smoothing the highest-frequency variations in the vertical water isotope profile. This diffusion process stops at the BCO depth at which water vapor can no longer move through pore space. The total amount of diffusion that occurs in the firn column depends on both the time it takes a parcel of firn to advect from the surface to the BCO depth and the temperature of the firn during that time. In analyses of water isotope records from ice cores, understanding this process of diffusion allows for both the correction of high-resolution water isotope records and interpretation of past firn conditions recorded by the ice core (Gkinis et al.2014; Jones et al.2017).

The CFM isotope diffusion module uses the equations presented in Johnsen et al. (2000). Each layer is assigned a water isotope δ value as a surface boundary condition. At each time step, the water isotope profiles δ(z,t) diffuse according to Fick's second law:

(26) D δ D t = Ω ( t , z ) 2 δ z 2 ,

where DδDt is the material derivative, and the isotope diffusivity Ω(t,z) depends on the firn temperature and density. This diffusive process smooths the profile; the profile is also affected by firn densification, which squeezes the δ values in adjacent layers together.

Both δ18O and δD can be tracked, accounting for differences in fractionation factors and air diffusivities. This module can be used to study cumulative water isotope diffusion under a range of firn conditions and/or to simulate water isotope records to compare to deep ice-core records. The CFM also tracks diffusion length.

2.5 Firn-air diffusion

The CFM includes a firn-air diffusion module coupled to the firn densification model. Previous modeling work has considered firn densification and firn-air transport separately; that is, firn-air models have assumed steady-state firn depth–density and effective diffusivity profiles. The firn-air module allows us to simulate gas transport while simultaneously modeling the evolution of firn depth and density in a changing climate.

The CFM firn-air module solves the firn-air equation (Severinghaus et al.2010; Birner et al.2018):

(27) C t = 1 ϕ op z ϕ op ( t , z ) κ eff ( t , z ) C z - Δ m g R T + Ω T z - w air ( z ) C z ,

where C is the concentration (ppm or ppt) or delta value (‰) of a gas species. The unitless parameter ϕop(t,z) is the open porosity; Δm is the molar mass difference between two isotopologues or the molar mass difference from air (kg mol−1); Ω is the thermal diffusion sensitivity (K−1), which is specific to individual gases (Severinghaus et al.2001); wair(z) is the advection rate of the air relative to the pore matrix; and κeff(t,z) is the effective diffusivity. As firn's porosity decreases, molecules must take a longer and more tortuous path as they diffuse. The effective diffusivity accounts for this by scaling the free-air diffusivity, κFA, with the firn's tortuosity τ (Buizert et al.2012): κeff=κFA/τ. In this work we take the term diffusivity to mean effective diffusivity.

Equation (27) is solved on the CFM's Lagrangian grid, which causes the downward advection of gases to be handled by the downward-moving reference frame. However, air advects downward slower than the surrounding firn. This occurs because the densification of the firn increases the air pressure in the open porosity, creating a pressure gradient. The air pressure gradient is positive downward, causing wair(z) to have a negative sign in a Lagrangian framework. Previous Lagrangian firn-air models have ignored this effect (e.g., Trudinger et al.1997), but it has been included in firn-air models with Eulerian grids (Buizert et al.2012); the CFM has the option to ignore it or to include it through the description provided by Buizert (2011).

The CFM also allows the user to choose which parameterization for diffusivity to use. As of the time of writing, the CFM includes the diffusivity parameterizations published by Schwander et al. (1988), Battle et al. (1996), Severinghaus et al. (2001), Freitag et al. (2002), Witrant et al. (2012), and Adolph and Albert (2014).

Modeling the diffusivity and the close-off physics requires knowing the BCO depth. To find the corresponding close-off density ρco, the CFM uses the relationship published by Martinerie et al. (1994):

(28) ρ co = 1 ρ ice + 6.95 × 10 - 7 T m - 4.3 × 10 - 5 - 1 ,

with temperature in Kelvin (K) and ρ in kilograms per cubic meter (kg m−3). Bubbles close off over a range of densities (and therefore depths). For example, Schwander and Stauffer (1984) reported that at Siple Station, 80 % of bubbles close off between 795 and 830 kg m−3. Equation (28) predicts the mean close-off density (and thus mean close-off porosity ϕco). The density of full bubble close-off, ρCO, where ϕop=0, will be slightly greater than ρco and its depth zco slightly deeper than z(ρco).

Equation (27) also requires knowing the open porosity ϕop, which equals the total porosity ϕtot minus the closed porosity ϕcl. The CFM uses the parameterization for ϕcl as a function of ϕtot presented in Goujon et al. (2003):

(29) ϕ cl = 0.37 ϕ tot ϕ tot ϕ co - 7.6 .

Finally, the CFM's firn-air module has an option to specify the lock-in depth (LID), the depth at which gravitational fractionation ceases despite the existence of open porosity. The lock-in zone (LIZ) is the zone between the LID and the close-off depth. The LIZ is not well-understood, but recent work has suggested that it may be created by the firn's three-dimensional layer structure and by barometric pumping (Birner et al.2018). In the CFM, the diffusivity below the specified LID is set to zero to inhibit further gas diffusion, and the lock-in density is determined by subtracting 14 kg m−3 from ρco (Blunier and Schwander2000).

The CFM does not include certain features that some firn-air models include (Buizert et al.2012) such as bubble trapping rate, bubble pressure, total air content, dispersive mixing in the LIZ, and the mean and distribution of gas ages in the closed porosity. These features will be integrated into future releases of the CFM.

2.6 Melt

The CFM has several meltwater percolation schemes to choose from, including two “tipping-bucket schemes”, a Richards equation single-domain scheme, and a Richards equation dual-domain approach. The latter two and one of the bucket schemes are described in Verjans et al. (2019). The second bucket scheme is similar to other bucket schemes that have been developed: at each time step, the volume of surface meltwater is allowed to percolate downward through the pore space. As the water reaches each model node (i.e., parcel of firn) in the model grid, the CFM first calculates the volume of water that refreezes due to the firn's cold content (the energy required to bring the firn's temperature to Tmelt), and that volume is immediately refrozen. The temperature of that parcel becomes the freezing temperature. Then, the volume of liquid that stays in the parcel due to capillary action (Coléou and Lesaffre1998) is subtracted from the meltwater volume. The remaining liquid moves downward to the next volume; this process continues until the entire volume of meltwater is accounted for. In the event that the firn's cold content can freeze the entire volume of meltwater, the firn temperature is raised by an amount that corresponds to the latent heat released by refreezing. If the meltwater encounters an impermeable layer, which we define as a layer with a density  800 kg m−3 (Gregory et al.2014), the water fills in the pore space in the parcel(s) above and remains liquid. After this percolation routine, the CFM solves for temperature in the entire firn column using the enthalpy scheme described in Sect. 2.3 and calculates the new mass and density of each parcel.

2.7 Grain growth and microstructure evolution

The CFM can optionally simulate grain size during a model run using one of two parameterizations, both of which assume spherical grains. The first gives the change in mean grain radius r (m) as in Gow (1969):

(30) d r d t = 1 2 r k g exp ( - E g / R T ) ,

with grain-growth activation energy Eg=42.4kJmol-1 and constant kg=1.3×10-7m2 s−1 taken from Cuffey and Paterson (2010, p. 40) and Arthern et al. (2010), respectively. Equation (30) is the grain-growth equation used for ART-T.

The second grain-growth parameterization accounts for the effect of liquid water on grain metamorphism (Sect. 2.6;  Verjans et al.2019). It is taken from Katsushima et al. (2009), who used equations from Tusima (1978) and Brun (1989) to simulate water flow in seasonal snowpacks. In the CFM, it is implemented as

(31) d r d t = 1 8 × 10 9 r 2 × min 2 π 1.28 × 10 - 8 + 4.22 × 10 - 10 θ 3 , 6.94 × 10 - 8 ,

where θ is mass-percent liquid water content.

As previously mentioned, the surface grain size of polar firn across the range of climates on ice sheets is not well-constrained. For Eq. (30), the CFM uses a uniform surface grain size r0=0.1mm. For Eq. (31), the CFM uses the empirical formula for surface grain size as a function of Tm (C) and b˙ (m w.e. a−1) given by Linow et al. (2012):

(32) r 0 = b 0 + b 1 T m + b 2 b ˙ ,

with constants b0=0.781, b1=0.0085, and b2=-2.79.

We recognize that many macroscale processes in firn (e.g., bulk densification) are dependent on the firn's microstructure (e.g., grain shape, size, and coordination number; specific surface area). Unfortunately, at present there is a lack of research describing the evolution of polar firn microstructure and how microstructure relates to macroscale firn processes. We have strived to design the CFM so that equations describing microscale evolution can be easily integrated into its framework. We will incorporate these equations when future research provides insights into how those properties evolve.

CFM applications

We demonstrate the utility of the CFM in two model applications. In the first, we compare the outputs of 13 firn densification models when forced with accumulation rates and temperatures predicted by a regional climate model (RCM). In the second, we use the coupled firn-density and firn-air model to simulate concentrations of gas stable isotopes trapped in ice cores during rapid climate changes.

3 Model application 1: intercomparison of firn model outputs at Summit, Greenland

In our first model application, we investigate uncertainty in firn model outputs that results from the choice of firn densification model. This work follows the Firn Model Intercomparison Experiment (FirnMICE; Lundin et al.2017), which compared the responses of eight firn densification models to synthetic climate histories that featured step changes in temperature and accumulation rate. Here, we use the CFM to expand upon that work by comparing outputs from 13 different firn densification models forced with temperature and accumulation histories for Summit, Greenland (72.58 N, 38.48 W; 3200 m). The mean annual temperature at Summit is −31.4C, and the annual accumulation averages 0.23 m ice eq. a−1. Historically, Summit has rarely experienced melt. Summit was the site of the GISP2 ice cores.

The FirnMICE project featured results submitted by different research groups running their own firn densification model codes; here, we run the firn densification models within the CFM framework. This allows us to compare the outputs from different firn densification models without concern about artifacts associated with different numerical methods (e.g., grid size, temperature diffusion scheme); that is, differences in model outputs are due to differences in the particular representation of physics in each firn densification model.

Sources of uncertainty in firn model outputs include the surface boundary conditions (i.e., the forcing) and the representation of physical processes in the firn densification model (i.e., the algorithms). We can begin to understand these uncertainties by comparing different outputs produced by a single firn densification model forced with range of plausible inputs or by comparing outputs from different firn densification models when they are forced by the same initial and boundary conditions. In this application we use the latter approach to leverage the CFM's ability to run multiple firn densification models. In particular, we examine the variability in model outputs that arises from firn densification model choice using the metrics of depth-integrated porosity (DIP), surface elevation change (dH), and bubble close-off (BCO) age and depth. In order to avoid picking a “best” firn densification model, which would only be best at our single test site, we focus on comparing the models' outputs to each other.

3.1 Methods

We forced each of the firn densification models in the CFM (Table 1) with skin temperature and surface mass balance (SMB) outputs for Summit, Greenland. from the RCM MAR3.9 (Fettweis et al.2017). We used the MAR products derived from ERA-40/ERA-Interim, which begin in January 1958 and end in October 2018. We ran the CFM at monthly time steps. The domain extended from the surface to  220 m of depth; the exact depth varied by model because the CFM's Lagrangian framework uses a fixed number of model nodes rather than a fixed grid. The predicted densities at that depth also varied by model but were generally 916 to 917 kg m−3. The surface density was held constant at 300 kg m−3, which is the mean value of the top 30 cm of 168 snow pits and firn cores from the Summit vicinity in the SUMup database (Montgomery et al.2018). For these simulations, we did not use the melt module and thus ignored any melt that may have occurred in the forcing fields during the simulation period.

In order to run a firn densification model, the model must first be spun up to an appropriate initial condition. Here, the initial condition that we desire is a firn depth–density–temperature profile for the start of the year 1958. Ideally, the spin-up process would produce a firn-density and firn-temperature profile that was the actual value of the profile in 1958 (i.e., what would have been measured in the field). At Summit, and most sites, that is not possible, so the goal is to create an initial condition that is as representative as possible of the firn at that time. To do this, we generated temperature and accumulation-rate histories for the 1000 years prior to the start of the model run (years 958 to 1957). Similar to Kuipers Munneke et al. (2015), we assumed that the 1958 to 1978 climate was in steady state and is representative of the climate for the previous 1000 years; the spin-up climate fields were made by repeating the 1958 to 1978 MAR temperature and SMB fields. We spun up the model for 1000 years because that was long enough to refresh the entire firn column (i.e., to remove any artifacts of the model initialization) and to ensure that the firn had come to thermal equilibrium. The firn during spin-up does not reach true steady state because there is climatological variability in the spin-up forcing data, but the firn does reach a state in which the variability in its properties (e.g., porosity) is consistent with a steady-state climate that includes natural variability.

3.2 Model intercomparison metrics

We compare the model results using several metrics. The depth-integrated porosity, DIP(z), is the volume of air contained within a 1×1m2 firn column (m3 m−2) above depth z, given by

(33) DIP ( z ) = 0 z ϕ ( z ) d z = 0 z ρ i - ρ ( z ) ρ i d z ,

where ϕ is the porosity, ρ(z) is the density at z, and ρi is the density of ice, taken in this work to be 917 kg m−3. DIP change is a key parameter used to convert volume change measurements (e.g., surface elevation change from altimetry) into mass change for sea-level rise estimates (Ligtenberg et al.2014). DIP is also called firn-air content (FAC).

When reporting the DIP predicted by a model, it is important to also report the maximum depth (and the corresponding density) to which the firn was modeled because if the bottom of the model domain is shallower than the transition to ice, there will be additional porosity beyond the model domain. In this study, we compare the models' predicted DIP in the upper 15 and 80 m and the total DIP (i.e., in the entire modeled firn column), which we refer to as DIP15, DIP80, and DIPtot, respectively.

When comparing the outputs from the firn densification models, one can consider both the inter-model differences in the predicted DIP and how the DIP changes through time in each model. The change in DIP is the quantity of interest for mass balance studies that need to adjust surface elevation measurements for interannual variability in firn thickness. Uncertainty in the total DIP is less consequential for these studies, but the total DIP is of interest when calculating the absolute mass of the ice sheet and predicting the volume of meltwater that can be retained on the ice sheet (e.g.,  Vandecrux et al.2019).

Predicting surface elevation change through time (dH∕dt) due to firn processes is essential for making corrections to surface elevation measurements from altimetry to derive mass changes. We calculate dH at each model time step by summing the ice-equivalent snow accumulation rate b˙, firn compaction rate vfc, vertical ice velocity vice at the bottom of the firn column due to dynamic ice-sheet processes, and vertical bedrock motion rate vbed:

(34) d H d t = b ˙ + v fc + v ice + v bed .

In steady state, the new snow accumulation rate is equal to the combined firn compaction and ice-sheet-thinning rates. In this application, we assume that vice is equal to the 1958–1978 mean ice-equivalent accumulation rate, which is in essence an assumption that the deep firn and ice sheet below are in steady state. We also assume for this application that vbed is zero and that additional layer thinning due to horizontal strain is negligible.

The BCO depth is the depth at which all porosity becomes closed (versus open and interconnected) and air is occluded in bubbles; the BCO age is the corresponding age of the firn at the BCO depth. The BCO density is commonly taken to be  830 kg m−3, but in reality, the BCO depth does not correspond to an exact density. Nevertheless, here we use the 830 kg m−3 density horizon for model comparisons. The BCO age is of interest to the ice-core science community as BCO age is a key parameter for determining Δage.

For each of these metrics, we calculate the mean, the standard deviation σ, and coefficient of variation (CV; the ratio of σ to the mean) of the models' results. There is no reason to believe that the mean of the models should give a better result (i.e., closer to observations) than a particular model; however, these statistics are useful for understanding how the models compare to one another.

Figure 1Profiles of depth–density (a, c) and depth–DIP (b, d) predicted by the models listed in Table 3 in the upper 80 m (a, b) and 15 m (c, d) of firn. Each panel is also included in the Supplement as a layered PDF file.


3.3 Model comparison results

3.3.1 Depth–density and depth–DIP profiles

Figure 1 shows depth–density (panels a, c) and depth–DIP(z) (panels b, d) profiles predicted by the various models in October 2018, the end of the model runs. Panels (a–b) and (c–d) show the model results to 80 and 15 m of depth, respectively. Each of the panels in Fig. 1 is included as a layered PDF file in the Supplement; the results from each model are plotted as a layer that can be toggled on and off using Adobe Reader or Acrobat software. Table 3 lists the models' predicted values of DIP15, DIP80, and DIPtot; the depth and age of the 830 kg m−3 density horizon (columns DEP830 and AGE830, respectively); and the linear trend of surface elevation change for the last 10 years of the model run. Table 3 also shows the arithmetic mean (“Model mean” row) and standard deviation (“Model σ” row) of the model results. These rows have two values for several columns; the values in parentheses are the statistics excluding CRO, which predicts anomalously low density at greater depths.

Because the depth and density of the firn are different at the bottom of the domain for each firn densification model, comparing the various models' DIPtot is not a direct comparison metric (as opposed to comparing DIP15). However, the amount of porosity that is in the bottom of the firn column is a very small percent of the total. For example, BAR reaches 916 kg m−3 at 126 m and 917 kg m−3 at 180 m; the DIP in this interval is 0.07 m, or 0.3 % of DIPtot. We thus consider DIPtot to be a worthy metric for comparison.

Near the surface (zone 1 densification), the models show a variety of responses. The mean DIP15 is 7.75 m with σ=0.81m, which gives CV =10.5 %. The DIP at 15 m of depth ranges from 6.4 m (LZ15, GOU) to 8.9 m (HEL). It is interesting that these three models at the high and low ends were all developed for Antarctica. If we exclude their results, the mean DIP15 becomes 7.89 m with σ=0.57 m. However, it is important to note that many sites in West Antarctica have a climate similar to that at Summit. We expect the models to show similar results near the surface because they all start with the same surface boundary condition (i.e., DIP(z=0m)=0m and ρsurface=300kg m−3). Additionally, one may expect models to perform well near the surface in the dry-firn zone because there is a relative abundance of shallow firn cores from Greenland and Antarctica that can be used for model calibration (see e.g., Kuipers Munneke et al.2015).

The depth–density results in zone 1 show how the models differ in their sensitivity to temperature. MW is the most sensitive to temperature, and it has the greatest density variability in zone 1. The GOU densification-rate equation for zone 1 is not a function of temperature, and CRO has relatively low temperature sensitivity. As such, these two models have very smooth depth–density profiles.

The models diverge in their predictions through zone 2. HEL and MW predict the highest DIP80, and ART-S predicts the lowest DIP80. The mean of the models' DIP80 is 20.8 m and σ is 2.99 m, giving a CV of 14.3 %.

Beyond 80 m, the models spread slightly more: excluding CRO, the CV is 16 %. The models should not be expected to spread significantly beyond 80 m because there is relatively little porosity beyond that depth, and all of the models are formulated to prevent densification beyond the ice density; i.e., they effectively have a fixed density boundary condition at the bottom of the domain.

CRO is a notable outlier in the deep firn. It predicts densification that is significantly slower than the other models. Its DIP80 is similar to the other models, but its density is ∼100kg m−3 less than the other models. CRO's DIPtot is 40 m, nearly twice the mean of the other models. This behavior is not surprising because CRO was developed as a seasonal snow model and is therefore not calibrated to accurately predict densification of higher-density firn. This is consistent with results from Lundin et al. (2017), who also found that a snow model did not predict densification well in deeper (zone 2) firn.

The mean depth of the 830 kg m−3 density horizon (excluding CRO) is 67 m, and σ=9.07m (CV =13 %). The mean age is 224 years with σ=29 years (CV =13 %). In general, models that predict deeper DEP830 also predict older AGE830 when compared with the others, but this is not universally true. For example, LZ15 predicts DEP830=65.46m and AGE830=232 years. LIG predicts a similar DEP830 of 64.37 m but a much younger AGE830 of 214 years (1 m of firn at Summit at this density is roughly 4 years of accumulation). Likewise, BAR predicts DEP830=70.03m and AGE830=231 years, while HEL predicts a deeper DEP830 (73.09 m) and younger AGE830 (227 years). These mismatches of BCO depths and ages demonstrate the effects of differences in model tuning. For example, consider a model that underestimates temperature sensitivity but overestimates b˙ sensitivity. In this case, these faulty sensitivities compensate for each other; the model will predict dρ∕dt and AGE830 reasonably well. However, due to the overestimated b˙ sensitivity, not enough material has been accumulated over the time required to reach a density of 830 kg m−3, and the model will predict a DEP830 that is too shallow.

Table 3Model results, including DIP at 15 and 80 m of depth and the bottom of the model domain ( 220–230 m); the depth and age of the 830 kg m−3 density horizon; and the linear trend in surface elevation change (i.e., a regression of the model results shown in Fig. 2) in the last 10 years of the model run (2008 to 2018). The “Model mean” and “Model σ” rows show the statistics for all 13 models; the values in parentheses are the statistics excluding CRO, which predicts anomalously low densification rates in the deeper firn.

Download Print Version | Download XLSX

3.3.2 Surface elevation change through time (dH∕dt)

Figure 2 shows the modeled surface elevation through time predicted by the models for the last 10 years of the model runs (October 2008 to October 2018). It is also included as a layered PDF file in the Supplement. Table 3 lists the linear least-squares trends (cm a−1) for each model from October in 2008–2016, 2016–2018, and 2008–2016 (columns dH∕dt08–16, dH∕dt16–18, and dH∕dt08–18, respectively). All the models predict that surface elevation has increased since 2008. The mean change from 2008 to 2018 is +12.5 cm, and the σ=2.7 cm (CV = 21 %).

In general, the models that predict the largest surface elevation increase consistently predict the largest increase throughout the entire time series and vice versa; i.e., the lines in Fig. 2 generally remain in the same location relative to one another. However, there are times that certain models change their relative positions. For example, the surface elevation increase since 2008 predicted by ART-S is the largest among the models as of mid-2016, but it is in the middle of the models at the end of the simulation. This behavior is a reflection of the models' different sensitivities to temperature, accumulation rate, and density.

Between 2008 and 2016, the change in surface elevation predicted by the models varies through time but does not deviate significantly from zero. The mean of the trends is −0.042cm a−1, and the models do not agree on the sign of the trend: six predict a small negative trend, and seven predict a small positive trend.

From October 2016 to the end of the model runs, the models all predict a surface elevation increase; this occurs because there are numerous months in that time period when MAR predicts that accumulation was higher than average. LZ15 predicts the largest trend (4.13cm a−1) and MW the smallest (2.61 cm a−1). The mean of the models' trends for these two years is 3.22 cm a−1 with σ=0.48cm a−1 (CV =15 %). The four models with the largest trend over this period (LZ15, HEL, HL, and LIG) were all tuned using Antarctic firn cores. There is a trade-off between sensitivity to temperature and sensitivity to accumulation rate involved in tuning a firn densification model, and the larger trends predicted by these “Antarctic” firn models may indicate that models tuned specifically for Antarctica are biased towards sensitivity to accumulation rate. For example, between LIG and KM, which are twin models tuned for Antarctica and Greenland, respectively, LIG predicts larger densification rates for the same mean accumulation rate (Eqs. 19 and 20) at sites with accumulation less than 0.8 m ice eq. a−1. This bias could occur if a large portion of the Antarctic cores came from sites with similar temperatures. Alternatively, models that are tuned for Greenland could be biased towards temperature sensitivity; e.g., MW, with the smallest trend over these 2 years, includes a significantly higher activation energy in the Arrhenius term.

Although the elevation change trend is clearly not linear over the October 2008 to 2018 period, fitting a linear trend to the modeled elevation changes further illustrates the differences between the models. In this case, the mean trend is 0.83 cm a−1, and σ=0.27cm a−1 (CV =32 %). ART-T predicts the smallest trend in dH∕dt (0.28cm a−1; the CV drops to 24 % if ART-T is excluded). MW predicts the smallest 2016–2018 trend, but it predicts the fourth-largest trend for 2008–2018. KM predicts the largest 2008–2018 trend (1.19cm a−1), whereas it had only the fifth-largest trend for 2016–2018.

Collectively, these results indicate that though the magnitude of surface elevation changes is relatively small, the models do not agree well with one another when simulating firn evolution in response to climate variability. The models predict different surface elevation trends relative to one another depending on the period considered, and in periods of relatively small changes in surface elevation, fitting trends to the modeled elevation change yields different signs depending on the model chosen.

3.3.3 Firn model uncertainty

Our results show that the firn densification model choice can be a significant source of uncertainty in applications requiring a firn model. Our goal with this application was to demonstrate the utility of the CFM in a simple model comparison exercise; as such, we have avoided detailed comparisons to data and instead focus on the broad agreement among the models. We ran 13 models, and they do not agree within 10 % when considering the DIP, BCO age and depth, or trend in surface elevation change. Models that agree well using one metric do not necessarily agree with a different metric. For example, KM and CRO predict nearly the same DIP15, but KM's 10-year trend in dH∕dt is the highest of all models, and CRO's is near the low end.

This is a challenge that the firn modeling community continues to face: despite the number of firn densification models that have been proposed, no single model is widely accepted. The general form of the firn densification models is relatively similar (e.g., dρ∕dt is a function of [ρiρ], which is an obvious “shut-off” to prevent over-densification), but they differ in their particular details (e.g., what the activation energy in the Arrhenius term is). Lundin et al. (2017) showed that firn densification models do not agree when predicting steady-state or transient behavior when forced with synthetic climate, and our results corroborate those results.

Future work could include an analysis of uncertainty related to firn model choice by running the suite of models across the entire range of the climates encountered on the ice sheets. There are additional sources of uncertainty in firn model outputs beyond the model choice; for example, any full uncertainty analysis will also require consideration of uncertainties in the boundary conditions and how those propagate through the model. The CFM is well-suited for such an exercise.

This model application focuses on a single location; the fact that the models do not agree well at Summit does not necessarily mean that they would not agree at other sites, but agreement is unlikely. If one model performs best at Summit, it does not necessarily indicate that model is the “best firn model”. Indeed, it may be the best model for Summit, but it may not work as well at other locations, and it is not obvious where one model might be “better” than another. A number of arguments could be made as to why it is inappropriate to compare these models, especially at a lone site in Greenland. For example, some of the models are potentially outdated and not in use any longer, certain models were tuned for a particular place and may not be appropriately applied to Summit, and some models were intended for ice-core Δage reconstruction rather than mass balance corrections or vice versa. Regardless, each of these models was at one point the state of the art, and each was designed to simulate the same properties of the firn. If the model equations are a correct and complete representation of the physics governing firn evolution, a model should be able to simulate firn evolution accurately on all timescales and all spatial scales.

Our results demonstrate a need to improve our understanding of firn densification physics, which may include both validation of existing models and development of new models. Unfortunately, data that are needed for the development of a purely physically based model are still lacking; any model development in the near future will require a certain amount of empirical tuning. For example, a microstructure-based firn densification model will need empirical parameterizations for the evolution of the microstructural properties. The addition of descriptions of physical processes such as grain growth to a model does not necessarily result in a better model if those physics (and the initial and boundary conditions) are not well-constrained. For example, ART-T includes grain growth, but it does not necessarily produce better results. Ultimately, research should be done to both (1) further our understanding of the microstructural evolution and underlying physics of firn evolution and (2) improve empirical models with observations of the firn's macroscale behavior.

Figure 2Model-predicted change in surface elevation at Summit, Greenland, from 2008 to 2018. This figure is also included in the Supplement as a layered PDF file.


4 Model application 2: firn-air stable isotopes and firn-thickness change during climate changes

In this application, we used the coupled firn-air and firn densification models to investigate the effect that a thickening and thinning (i.e., non-steady-state) firn column has on nitrogen and argon isotope records in ice cores.

Figure 3(a) Accumulation-rate and temperature histories derived from the GISP2 ice core (Alley2004; Cuffey and Clow1997) and used to force the CFM for the analyses in Sect. 4. (b) LID predicted by the steady-state and transient model simulations described in Sect. 4. The LID for the transient simulation is determined by running the CFM with the climate shown in (a), and the LID for the steady-state simulation is that predicted by the Herron and Langway (1980) model with the accumulation rate 0.07 m ice eq. a−1 and mean annual temperature −47.5C.


Gas isotopes fractionate in firn due to gravity (heavier isotopes become enriched at greater depths) and due to thermal gradients (heavier isotopes become enriched in the colder part of the firn column) (Severinghaus et al.1998). In a steady climate, there are seasonal temperature gradients near the surface (upper  10 m of firn), and the deeper firn is isothermal or near isothermal due to the low temperature diffusivity of firn. During climate change events, the surface warms or cools and creates a temperature gradient between the surface and the lock-in depth (LID), causing thermal fractionation. Different gas species have different thermal sensitivities, and this fact can be leveraged to infer the magnitude of temperature changes during rapid climate change events. Previous work showed that the temperature at Summit increased by  9 C over the course of several decades during the Bølling transition (14.67 ka before 1950; Severinghaus and Brook1999) and by 5–10 C over a century at the end of the Younger Dryas (11.6 ka before 1950;  Severinghaus et al.1998). In both these cases, the authors examined δ15N and δ40Ar isotopes; δ15N and δ40Ar∕4 will have the same gravitational fractionation signal. Any deviation of these species from one another in the firn is a result of thermal fractionation. Severinghaus and Brook (1999) and Severinghaus et al. (1998) modeled gas isotopes to infer temperature increases, and part of their data–model mismatch was attributed to the fact that their model assumed a steady-state firn column. This model was unable to account for transient firn thickening due to an accumulation-rate increase coincident with the temperature increase; their model predicted values  3 ‰–4 ‰ less than the observed values.

4.1 Methods

We used the CFM's coupled firn-air and firn-density models to examine the effect of transient firn evolution on gas records in ice cores during rapid climate change events. We ran two model simulations of the evolution of δ15N and δ40Ar in the firn at Summit. We forced the CFM with temperature and accumulation-rate histories from the GISP2 ice core (Cuffey and Clow1997; Alley2000, 2004). Both simulations ran for 49 000 years, which is the length of the climate records. We ran the model using yearly time steps, and the model domain extended to a depth of  2200 m. The difference between the two model runs was the firn depth–density profile: in the first simulation, we used a constant profile; in the second, we used a transient firn densification model to allow the density to evolve with the climate.

For the first simulation, we used a steady-state depth–density profile predicted by the Herron and Langway (1980) analytic model with an accumulation rate of 0.07 m ice eq. a−1 and a temperature of −47.5C, which are consistent with values during the Younger Dryas and leading into the Bølling transition. They are also the values used for modeling in Severinghaus et al. (1998) and Severinghaus and Brook (1999). In this simulation, we used the GISP2 temperature record (Fig. 3a) to force the CFM's temperature evolution module, thereby allowing the firn temperature to evolve. However, we did not allow this temperature forcing to affect the model depth–density profile, and the LID stayed constant (Fig. 3b) at  97 m.

For the second simulation, we used the GISP2 temperature and accumulation-rate data (Fig. 3, upper panel) to run the CFM in transient mode using the Herron and Langway (1980) dynamic firn densification model. In addition to the firn-temperature profile evolution, this simulation allowed the firn depth–density profile to evolve, causing the LID to vary through time (Fig. 3b).

Figure 4Measured and modeled δ15N and δ40Ar∕4 profiles. Data are from Severinghaus et al. (1998) and Severinghaus and Brook (1999). We add 1.5 % to the modeled gas ages to fit the data.


4.2 Firn-air results

Figure 4 shows the results of the two simulations and the δ15N and δ40Ar∕4 data from Severinghaus et al. (1998) and Severinghaus and Brook (1999). The horizontal axis of the plot is the gas age. We add 1.5 % to the modeled gas ages because the gas ages predicted by the model are too young compared to the data (determined by comparing the timing of the modeled isotope increases to the data during the Bølling transition and Younger Dryas). The model is likely failing to produce the correct gas age for several reasons: (1) for our simple experiment, we assumed a uniform gas age of 15 years at the LID; (2) the large modeled isotope changes occur when the firn densification model, which is not necessarily accurate, predicts a firn-thickness change; and (3) there are likely uncertainties in the climate forcing data, including the timescale of those data. The δ15N and δ40Ar∕4 data were converted from depth to gas age using the GICC05 timescale (Rasmussen et al.2014; Seierstad et al.2014).

For the steady-state simulations, variability in the isotope values is due only to fractionation from temperature gradients in the firn. For the transient model runs, the variability is due to both fractionation from temperature gradients and to changes in the firn-column thickness. For the isotopes considered independently, much of the difference in the predicted isotope values between the transient and steady-state simulations shown in Fig. 4 can be attributed to the change in firn-column thickness. For example, at  14 500 years before present, the peak of the δ15N predicted by the transient model is about 0.02 ‰ higher than the steady-state run. In this case, the increased firn thickness results in more gravitational fractionation. However, the transient and steady curves are also offset from one another temporally. This is because (1) the temperature gradients that form in the transient firn are different from the gradients in the steady-state firn; i.e., a thicker diffusive column in the firn will result in a different temperature gradient through time for the same surface temperature increase, and (2) the timescales of diffusion are slightly different; i.e., it will take longer for the thicker firn to come to a new thermal equilibrium.

Figure 5As in Fig. 4, but zoomed into the Bølling transition. Additionally, for this figure we have shifted the y axis to match the transient and steady-state simulations with the data leading into the Bølling transition to highlight the magnitude of the modeled changes compared to the data.


On the whole, the transient model matches the data better than the steady-state model, especially during the Younger Dryas and its termination. Notably, the transient model does well at predicting the high values of δ15N and δ40Ar∕4 associated with the rapid warming at the end of the Younger Dryas. These values are  0.03 ‰ higher than those predicted by the steady-state model and are consistent with the model–data misfit in Severinghaus et al. (1998), who used a fixed LID.

There is a  0.05 ‰ offset between the data and transient model at  14.8 ka, suggesting that the model is predicting an LID that is  10 m too shallow ( 10 % of the firn-column thickness) at the start of the Bølling transition. This may be caused by forcing data uncertainty: at around −47.5C and 0.07 m ice eq. a−1, the Herron and Langway (1980) model predicts a 4 m change in LID for a 1 C change in temperature and a 5 m change in LID for a 0.01 m ice eq. a−1 change in accumulation rate; a small uncertainty in the forcing data for either of these variables could account for the 0.05 ‰ offset. The model–data discrepancy could also be caused by model inadequacy: the Herron and Langway (1980) model's calibration data set included only two cold low-accumulation sites, so its accuracy in those conditions may be questionable. The model also may not be predicting a large enough temperature gradient, which could occur if the temperature increase in the forcing data was not large enough.

Figure 5 shows the data and model results zoomed in on the Bølling transition. In order to directly compare the magnitude of the modeled δ15N and δ40Ar∕4 to the data, we have shifted the data and the steady-state and transient isotope values by subtracting the mean 15–14.75 ka δ15N values. The transient model matches the observed magnitude and rate of the δ15N and δ40Ar∕4 change much better than the steady-state model. This is because the transient CFM predicts that the LID increases from 92 m at 14.7 ka to 105 m at 14.4 ka.

Despite the lingering model–data mismatch, the transient model clearly performs better than the steady-state model, and our model results during the Younger Dryas and Bølling transition corroborate the assertion in Severinghaus and Brook (1999) that their model–data misfit is due to transient firn thickening. Our results do not change their conclusions but do provide assurance that their conclusions are sound. This application also demonstrates the utility of the CFM's coupled firn-air and firn-density model. In this case, the coupled model should allow a more accurate assessment of the magnitude of temperatures because it can account for the temperature gradient that will form in thickening firn, which is smaller than that in steady-state firn. It is important to note that most firn densification models (including Herron and Langway1980) were developed with a steady-state assumption, and applying them to transient simulations produces additional uncertainty. This uncertainty is challenging to quantify, however, because we do not have direct observations of firn evolution during rapid climate changes.

The CFM's coupled firn densification and firn-air modules have additional potential to help test hypotheses surrounding anomalies in ice-core records. For example, experiments could be done to investigate the impact that an impermeable ice lens would have on gas records. It could also be used to model water isotope diffusion simultaneously with firn-air transport. The CFM's ability to model multiple physical processes in a single framework allows us to investigate processes with different timescales. For example, a temperature change (with no concurrent accumulation-rate change) will create a temperature gradient in the firn and will also cause the firn to change thickness by affecting the densification rate, but those processes will operate on different timescales.

4.3 Additional CFM applications

The CFM been used in several studies to date, and we briefly mention those here. Verjans et al. (2019) used the CFM to compare model outputs from three firn meltwater percolation schemes using HL, KM, and CRO. The firn meltwater Retention Model Intercomparison Project (RetMIP) (Vandecrux et al.2020) compared results from nine firn densification and hydrology models, including the CFM. Garland et al. (2018) used the CFM's isotope diffusion module to examine the correlation between accumulation and the weighted permutation entropy of isotopes in ice-core records. Hughes et al. (2020) used the isotope diffusion module to help interpret the water isotope records from an ice core drilled on the Renland ice cap, Greenland.

5 Conclusions

We developed the Community Firn Model (CFM), an open-source firn model framework. The CFM includes modules to simulate a number of physical processes in firn, including densification, heat transport, meltwater percolation, grain growth, water isotope diffusion, and firn-air diffusion. We demonstrated the utility of the CFM in two model applications. In the first, we leveraged the CFM's ability to run numerous firn densification models by forcing 13 models with the regional climate model outputs for Summit, Greenland. These simulations showed that the choice of firn densification model can contribute significant uncertainty to firn model outputs: the model spread was greater than 10 % of the model mean for the metrics of depth-integrated porosity, bubble close-off depth and age, and surface elevation change trend. There is no single densification model that is widely considered best; different models are preferred for different locations (e.g., Antarctica vs. Greenland) or for different applications (e.g., ice cores vs. satellite altimetry). Continued studies are necessary to improve firn densification models and to better understand the uncertainty in their applications. These include investigations of the microstructural evolution of firn and in situ measurements of the bulk firn densification rate in a variety of climates.

In the second application, we investigated the effect of a thickening or thinning firn column on noble gas isotope records in ice cores. To our knowledge, the CFM is the first model that couples transient firn densification and firn-air transport. This application demonstrated that the coupled model can predict records of isotopes in ice cores better than a firn-air model that uses a steady-state firn-density profile. This tool could be used for a number of studies surrounding firn-air transport in changing climates. The model is limited in that it relies on published parameterizations of the effective diffusivity, which may fail to incorporate all relevant parameters, especially near and in the lock-in zone. Continued research investigating microstructure at the bottom of the firn column is needed to improve our ability to model air transport accurately.

The goals of the CFM project are to provide a community resource that can be used by research groups that need a firn model and to provide the ability to create open-source results for model comparison and benchmarking. The CFM has already been used for several studies, including Verjans et al. (2019) and Garland et al. (2018). The CFM allows a fast and easy way to run a model experiment using the same boundary conditions with different densification physics, and it removes potential sources of “noise” when comparing the outputs of different models, e.g., numerical solvers and different time stepping. The code is open source, allowing anyone to be able to check results from researchers who use the CFM. As the firn research community improves our understanding of physical processes in firn, e.g., new descriptions of densification or meltwater percolation processes, the new knowledge can be incorporated easily into the CFM's modular framework. In the spirit of open-source software, we encourage other researchers to develop their own modules and add them to the code or to request features that would improve the CFM's usefulness and/or ease of use.

Code and data availability

The CFM code is publicly available under the MIT license at (last access: 10 September 2020 Stevens et al.2019). Its documentation is online at (last access: 10 September 2020). All model outputs and scripts used to make the figures are freely available upon request.


The supplement related to this article is available online at:

Author contributions

CMS was the main developer of the CFM, conceived and ran the model applications, and wrote the paper. EDW initiated the CFM project and supervised its development. JL helped initiate the CFM project, planned model architecture, and developed code. VV, EK, AH, and BH developed code for the CFM. All authors contributed to paper writing and editing.

Competing interests

The authors declare that they have no conflict of interest.


We thank Paul Harris, Alex Le, Will Leahy, Huong Vo, and Michael Yoon for providing valuable assistance developing the CFM's code. Thank you to Brooke Medley, Vasileios Gkinis, and Jiajen Chen for providing feedback and bug reports on the CFM. Thank you to Tyler J. Fudge, Michelle Koutnik, and Howard Conway for constructive conversations that improved the CFM and this paper. We used GNU Parallel (Tange2011) to parallelize model runs for application 1.

We also thank two anonymous reviewers for providing useful critiques that improved the paper.

Financial support

This research has been supported by the National Science Foundation (grant nos. 0968391 and 1443471) and the National Aeronautics and Space Administration (grant no. NNX15AC62G).

Review statement

This paper was edited by Philippe Huybrechts and reviewed by two anonymous referees.


Adolph, A. C. and Albert, M. R.: Gas diffusivity and permeability through the firn column at Summit, Greenland: measurements and comparison to microstructural properties, The Cryosphere, 8, 319–328,, 2014. a

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,, 2019. a

Alexander, P., Tedesco, M., Koenig, L., and Fettweis, X.: Evaluating a regional climate model simulation of Greenland ice sheet snow and firn density for improved surface mass balance estimates, Geophys. Res. Lett., 46, 12073–12082,, 2019. a

Alley, R. B.: Firn densification by grain-boundary sliding: a first model, Le Journal de Physique Colloques, 48, C1–249, 1987. a, b, c, d

Alley, R. B.: The Younger Dryas cold interval as viewed from central Greenland, Quatern. Sci. Rev., 19, 213–226,, 2000. a

Alley, R. B.: GISP2 Ice Core Temperature and Accumulation Data, IGBP PAGES/World Data Center for Paleoclimatology. Data Contribution Series #2004-013, NOAA/NGDC Paleoclimatology Program, Boulder CO, USA, 2004. a, b

Anderson, E. A.: A point energy and mass balance model of a snow cover, NOAA technical report NWS; 19, Office of Hydrology, National Weather Service, 1976. a

Arnaud, L., Barnola, J. M., and Duval, P.: Physical modeling of the densification of snow/firn and ice in the upper part of polar ice sheets, in: Physics of ice core records, Hokkaido University Press, 285–305, 2000. a, b, c

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, J. Geophys. Res.-Ea. Surf., 115, F03011,, 2010. a, b, c, d, e, f, g

Arzt, E.: The influence of an increasing particle coordination on the densification of spherical powders, Acta Metall. Mater., 30, 1883–1890, 1982. a, b

Bader, H.: Sorge's Law of Densification of Snow on High Polar Glaciers, J, Glaciol,, 2, 319–323,, 1954. a

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,, 1991. a, b, c

Battle, M., Bender, M., Sowers, T., Tans, P. P., Butler, J. H., Elkins, J. W., Ellis, J. T., Conway, T., Zhang, N., Lang, P., and Clarket, A. D.: Atmospheric gas concentrations over the past century measured in air from firn at the South Pole, Nature, 383, 231–235,, 1996. a

Birner, B., Buizert, C., Wagner, T. J. W., and Severinghaus, J. P.: The influence of layering and barometric pumping on firn air transport in a 2-D model, The Cryosphere, 12, 2021–2037,, 2018. a, b

Blunier, T. and Schwander, J.: Gas enclosure in ice: age difference and fractionation, in: Physics of Ice Core Records, Hokkaido University Press, 307–326, 2000. a, b

Brandt, R. E. and Warren, S. G.: Temperature measurements and heat transfer in near-surface snow at the South Pole, J. Glaciol., 43, 339–351,, 1997. a

Brun, E.: Investigation on Wet-Snow Metamorphism in Respect of Liquid-Water Content, Ann. Glaciol., 13, 22–26,, 1989. a

Brun, E., David, P., Sudul, M., and Brunot, G.: A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting, J. Glaciol., 38, 13–22,, 1992. a

Buizert, C.: The influence of firn air transport processes and radiocarbon production on gas records from polar firn and ice, PhD thesis, University of Copenhagen, 2011. a

Buizert, C., Martinerie, P., Petrenko, V. V., Severinghaus, J. P., Trudinger, C. M., Witrant, E., Rosen, J. L., Orsi, A. J., Rubino, M., Etheridge, D. M., Steele, L. P., Hogan, C., Laube, J. C., Sturges, W. T., Levchenko, V. A., Smith, A. M., Levin, I., Conway, T. J., Dlugokencky, E. J., Lang, P. M., Kawamura, K., Jenk, T. M., White, J. W. C., Sowers, T., Schwander, J., and Blunier, T.: Gas transport in firn: multiple-tracer characterisation and model intercomparison for NEEM, Northern Greenland, Atmos. Chem. Phys., 12, 4259–4277,, 2012. a, b, c, d, e

Buizert, C., Cuffey, K. M., Severinghaus, J. P., Baggenstos, D., Fudge, T. J., Steig, E. J., Markle, B. R., Winstrup, M., Rhodes, R. H., Brook, E. J., Sowers, T. A., Clow, G. D., Cheng, H., Edwards, R. L., Sigl, M., McConnell, J. R., and Taylor, K. C.: The WAIS Divide deep ice core WD2014 chronology – Part 1: Methane synchronization (68–31 ka BP) and the gas age–ice age difference, Clim. Past, 11, 153 173,, 2015. a, b, c

Coléou, C. and Lesaffre, B.: Irreducible water saturation in snow: experimental results in a cold laboratory, Ann. Glaciol., 26, 64–68,, 1998. a

Cuffey, K. M. and Clow, G. D.: Temperature, accumulation, and ice sheet elevation in central Greenland through the last deglacial transition, J. Geophys. Res.-Oceans, 102, 26383–26396,, 1997. a, b

Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, Elsevier, 2010. a

Cullather, R. I., Nowicki, S. M. J., Zhao, B., and Koenig, L. S.: A Characterization of Greenland Ice Sheet Surface Melt and Runoff in Contemporary Reanalyses and a Regional Climate Model, Front. Earth Sci., 4, 10 pp.,, 2016. a

Fausto, R. S., Box, J. E., Vandecrux, 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, Front. Earth Sci., 6, 51 pp.,, 2018. a

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,, 2017. a, b, c

Freitag, J., Dobrindt, U., and Kipfstuhl, J.: A new method for predicting transport properties of polar firn with respect to gases on the pore-space scale, Ann. Glaciol., 35, 538–544,, 2002. a

Garland, J., Jones, T. R., Bradley, E., Neuder, M., and White, J. W. C.: Climate entropy production recorded in a deep Antarctic ice core, 15 pp., arXiv:1806.10936, 2018. a, b

Gkinis, V., Simonsen, S., Buchardt, S., White, J., and Vinther, B.: Water isotope diffusion rates from the NorthGRIP ice core for the last 16,000 years – Glaciological and paleoclimatic implications, Earth Planet. Sci. Lett., 405, 132–141,, 2014. a

Goujon, C., Barnola, J.-M., and Ritz, C.: Modeling the densification of polar firn including heat diffusion: Application to close-off characteristics and gas isotopic fractionation for Antarctica and Greenland sites, J. Geophys. Res.-Atmos., 108, 4792,, 2003. a, b, c, d, e, f, g, h, i

Gow, A. J.: On the Rates of Growth of Grains and Crystals in South Polar Firn, J. Glaciol., 8, 241–252,, 1969. a

Gow, A. J.: Time-temperature dependence of sintering in perennial isothermal snowpacks, in: IAHS-AISH Publication No. 114, 25–41, 1975. a

Gow, A. J., Meese, D. A., and Bialas, R. W.: Accumulation variability, density profiles and crystal growth trends in ITASE firn and ice cores from West Antarctica, Ann. Glaciol., 39, 101–109,, 2004. a

Gregory, S. A., Albert, M. R., and Baker, I.: Impact of physical properties and accumulation rate on pore close-off in layered firn, The Cryosphere, 8, 91–105,, 2014. a

Helsen, M. M., van den Broeke, M. R., van de Wal, R. S. W., van de Berg, W. J., van Meijgaard, E., Davis, C. H., Li, Y., and Goodwin, I.: Elevation Changes in Antarctica Mainly Determined by Accumulation Variability, Science, 320, 1626–1629,, 2008. a, b, c

Herron, M. M. and Langway, C. C.: Firn Densification: An Empirical Model, J. Glaciol., 25, 373–385,, 1980. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q

Hughes, A. G., Jones, T. R., Vinther, B. M., Gkinis, V., Stevens, C. M., Morris, V., Vaughn, B. H., Holme, C., Markle, B. R., and White, J. W. C.: High-frequency climate variability in the Holocene from a coastal-dome ice core in east-central Greenland, Clim. Past, 16, 1369–1386,, 2020. a

Hulbe, C. L. and Whillans, I. M.: A method for determining ice-thickness change at remote locations using GPS, Annal. Glaciol., 20, 263–268,, 1994. a

Jiawen, R., Dahe, Q., and Maohuan, H.: Thermal properties and temperature distribution of snow/firn on the Law Dome ice cap, Antarctica, Chinese J. Polar Sci., 2, 38–46, 1991. a

Johnsen, S. J., Clausen, H. B., Cuffey, K. M., Hoffmann, G., Schwander, J., and Creyts, T.: Diffusion of stable isotopes in polar firn and ice: the isotope effect in firn diffusion, in: Physics of ice core records, Hokkaido University Press, 121–140, 2000. a

Jones, T., Cuffey, K., White, J., Steig, E., Buizert, C., Markle, B., McConnell, J., and Sigl, M.: Water isotope diffusion in the WAIS Divide ice core during the Holocene and last glacial, J. Geophys. Res.-Ea. Surf., 122, 290–309,, 2017. a

Katsushima, T., Kumakura, T., and Takeuchi, Y.: A multiple snow layer model including a parameterization of vertical water channel process in snowpack, Cold Reg. Sci. Technol., 59, 143–151,, 2009. a

Kawamura, K., Severinghaus, J. P., Ishidoya, S., Sugawara, S., Hashida, G., Motoyama, H., Fujii, Y., Aoki, S., and Nakazawa, T.: Convective mixing of air in firn at four polar sites, Earth Planet. Sc. Lett., 244, 672–682,, 2006. a

Kuipers Munneke, P., van den Broeke, M. R., Reijmer, C. H., Helsen, M. M., Boot, W., Schneebeli, M., and Steffen, K.: The role of radiation penetration in the energy budget of the snowpack at Summit, Greenland, The Cryosphere, 3, 155–165,, 2009. a

Kuipers Munneke, P., Ligtenberg, S. R. M., Noël, B. P. Y., Howat, I. M., Box, J. E., Mosley-Thompson, E., McConnell, J. R., Steffen, K., Harper, J. T., Das, S. B., and van den Broeke, M. R.: Elevation change of the Greenland Ice Sheet due to surface mass balance and firn processes, 1960–2014, The Cryosphere, 9, 2009–2025,, 2015. a, b, c, d, e, f, g, h

Langen, P. L., Fausto, R. S., Vandecrux, B., Mottram, R. H., and Box, J. E.: Liquid water flow and retention on the Greenland ice sheet in the regional climate model HIRHAM5: Local and large-scale impacts, Front. Earth Sci., 4, 110 pp.,, 2017. a, b

Li, J. and Zwally, H. J.: Modeled seasonal variations of firn density induced by steady-state surface air-temperature cycle, Ann. Glaciol., 34, 299–302,, 2002. a, b

Li, J. and Zwally, H. J.: Modeling the density variation in the shallow firn layer, Anna. Glaciol., 38, 309–313, 2004. a

Li, J. and Zwally, H. J.: Modeling of firn compaction for estimating ice-sheet mass change from observed ice-sheet elevation change, Ann. Glaciol., 52, 1–7,, 2011. a, b, c, d

Li, J. and Zwally, H. J.: Response times of ice-sheet surface heights to changes in the rate of Antarctic firn compaction caused by accumulation and temperature variations, J. Glaciol., 61, 1037–1047,, 2015. a, b, c

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,, 2011. a, b, c, d

Ligtenberg, S. R. M., Kuipers Munneke, P., and van den Broeke, M. R.: Present and future variations in Antarctic firn air content, The Cryosphere, 8, 1711–1723,, 2014. a

Linow, S., Hörhold, M., and Freitag, J.: Grain-size evolution of polar firn: A new empirical grain growth parameterization based on X-ray microcomputer tomography measurements, J. Glaciol., 58, 1245–1252,, 2012. a

Lundin, J. M., Stevens, C. M., Arthern, R., Buizert, C., Orsi, A., Ligtenberg, S. R., Simonsen, S. B., Cummings, E., Essery, R., Leahy, W., Harris, P., Helsen, M. M., and Waddington, E. D.: Firn Model Intercomparison Experiment (FirnMICE), J. Glaciol., 63, 401–422,, 2017. a, b, c, d, e, f

Lüthi, M. P. and Funk, M.: Modelling heat flow in a cold, high-altitude glacier: interpretation of measurements from Colle Gnifetti, Swiss Alps, J. Glaciol., 47, 314–324, 2001. a

Maeno, N. and Ebinuma, T.: Pressure sintering of ice and its implication to the densification of snow at polar glaciers and ice sheets, J. Phys. Chem., 87, 4103–4110,, 1983. a

Martinerie, P., Lipenkov, V. Y., Raynaud, D., Chappellaz, J., Barkov, N., and Lorius, C.: Air content paleo record in the Vostok ice core (Antarctica): A mixed record of climatic and glaciological parameters, J. Geophys. Res.-Atmos., 99, 10565–10576,, 1994. a

Meyer, C. R. and Hewitt, I. J.: A continuum model for meltwater flow through compacting snow, The Cryosphere, 11, 2799–2813,, 2017. a

Montgomery, L., Koenig, L., and Alexander, P.: The SUMup dataset: compiled measurements of surface mass balance components over ice sheets and sea ice with analysis over Greenland, Earth Syst. Sci. Data, 10, 1959–1985,, 2018. a

Morris, E. M. and Wingham, D. J.: Densification of polar snow: Measurements, modeling, and implications for altimetry, J. Geophys. Res.-Ea. Surf., 119, 349–365,, 2014. a, b, c, d, e, f, g, h

Noël, B., van de Berg, W. J., van Wessem, J. M., van Meijgaard, E., van As, D., Lenaerts, J. T. M., Lhermitte, S., Kuipers Munneke, P., Smeets, C. J. P. P., van Ulft, L. H., van de Wal, R. S. W., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 1: Greenland (1958–2016), The Cryosphere, 12, 811–831,, 2018. a, b

Patankar, S. V.: Numerical heat transfer and fluid flow, CRC Press, Boca Raton,, 1980. a

Pimienta, P.: Etude du comportement mécanique des glaces polycristallines aux faibles contraintes: applications aux glaces des calottes polaires, PhD thesis, Université Scientifique, Technique et Médicale de Grenoble, 1987. a

Rasmussen, S. O., Bigler, M., Blockley, S. P., Blunier, T., Buchardt, S. L., Clausen, H. B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S. J., Fischer, H., Gkinis, V., Guillevic, M., Hoek, W. Z., Lowe, J. J., Pedro, J. B., Popp, T., Seierstad, I. K., Steffensen, J. P., Svensson, A. M., Vallelonga, P., Vinther, B. M., Walker, M. J., Wheatley, J. J., and Winstrup, M.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy, Quatern. Sci. Rev., 106, 14–28,, 2014. a

Riche, F. and Schneebeli, M.: Thermal conductivity of snow measured by three independent methods and anisotropy considerations, The Cryosphere, 7, 217–227,, 2013. a

Robin, G. d. Q.: Glaciology III: Seismic shooting and related investigations, in: Norwegian-British-Swedish Antarctic Expedition, 1949–52, Scientific Results, vol. 5, Norsk Polarinstit, Oslo, 1958. a

Schwander, J. and Stauffer, B.: Age difference between polar ice and the air trapped in its bubbles, Nature, 311, 45–47,, 1984. a

Schwander, J., Stauffer, B., and Sigg, A.: Air Mixing in Firn and the Age of the Air at Pore Close-Off, Ann. Glaciol., 10, 141–145,, 1988. a

Schwander, J., Sowers, T., Barnola, J.-M., Blunier, T., Fuchs, A., and Malaizé, B.: Age scale of the air in the summit ice: Implication for glacial-interglacial temperature change, J. Geophys. Res.-Atmos., 102, 19483–19493,, 1997. a

Schwerdtfeger, P.: Theoretical derivation of the thermal conductivity and diffusivity of snow, Int. Assoc. Sci. Hydrol. Publ, 61, 75–81, 1963. a

Schytt, V.: Glaciology II: Snow Studies at Maudheim: Snow Studies Inland: the Inner Structure of the Ice Shelf at Maudheim as Shown by Core Drilling, in: Norwegian-British-Swedish Antarctic Expedition, 1949–52, Scientific Results, vol. 4, Norsk Polarinstitutt, Oslo, 1958. a

Seierstad, I. K., Abbott, P. M., Bigler, M., Blunier, T., Bourne, A. J., Brook, E., Buchardt, S. L., Buizert, C., Clausen, H. B., Cook, E., Dahl-Jensen, D., Davies, S. M., Guillevic, M., Johnsen, S. J., Pedersen, D. S., Popp, T. J., Rasmussen, S. O., Severinghaus, J. P., Svensson, A., and Vinther, B. M.: Consistently dated records from the Greenland GRIP, GISP2 and NGRIP ice cores for the past 104 ka reveal regional millennial-scale δ18O gradients with possible Heinrich event imprint, Quatern. Sci. Rev., 106, 29– 46,, 2014. a

Severinghaus, J. P. and Brook, E. J.: Abrupt Climate Change at the End of the Last Glacial Period Inferred from Trapped Air in Polar Ice, Science, 286, 930–934,, 1999. a, b, c, d, e, f

Severinghaus, J. P., Sowers, T., Brook, E. J., Alley, R. B., and Bender, M. L.: Timing of abrupt climate change at the end of the younger dryas interval from thermally fractionated gases in polar ice, Nature, 391, 141–146,, 1998. a, b, c, d, e, f, g, h

Severinghaus, J. P., Grachev, A., and Battle, M.: Thermal fractionation of air in polar firn by seasonal temperature gradients, Geochemistry, Geophysics, Geosystems, 2, 1048,, 2001. a, b

Severinghaus, J. P., Albert, M. R., Courville, Z. R., Fahnestock, M. A., Kawamura, K., Montzka, S. A., Mühle, J., Scambos, T. A., Shields, E., Shuman, C. A., Suwa, M., Tans, P., and Weiss, R. F.: Deep air convection in the firn at a zero-accumulation site, central Antarctica, Earth Planet. Sci. Lett., 293, 359–367,, 2010. a

Shepherd, A., Ivins, E. R., A, G., Barletta, V. R., Bentley, M. J., Bettadpur, S., Briggs, K. H., Bromwich, D. H., Forsberg, R., Galin, N., Horwath, M., Jacobs, S., Joughin, I., King, M. A., Lenaerts, J. T. M., Li, J., Ligtenberg, S. R. M., Luckman, A., Luthcke, S. B., McMillan, M., Meister, R., Milne, G., Mouginot, J., Muir, A., Nicolas, J. P., Paden, J., Payne, A. J., Pritchard, H., Rignot, E., Rott, H., Sørensen, L. S., Scambos, T. A., Scheuchl, B., Schrama, E. J. O., Smith, B., Sundal, A. V., van Angelen, J. H., van de Berg, W. J., van den Broeke, M. R., Vaughan, D. G., Velicogna, I., Wahr, J., Whitehouse, P. L., Wingham, D. J., Yi, D., Young, D., and Zwally, H. J.: A Reconciled Estimate of Ice-Sheet Mass Balance, Science, 338, 1183–1189,, 2012. a

Simonsen, S. B., Stenseng, L., Adalgeirsdottir, G., Fausto, R. S., Hvidberg, C. S., and Lucas-Picher, P.: Assessing a multilayered dynamic firn-compaction model for Greenland with ASIRAS radar measurements, J. Glaciol., 59, 545–558,, 2013. a, b, c, d

Sowers, T., Bender, M., Raynaud, D., and Korotkevich, Y. S.: δ15N of N2 in air trapped in polar ice: A tracer of gas transport in the firn and a possible constraint on ice age-gas age differences, J. Geophys. Res.-Atmos., 97, 15683–15697,, 1992. a

Steffen, K. and Box, J.: Surface climatology of the Greenland Ice Sheet: Greenland Climate Network 1995–1999, J. Geophys. Res.-Atmos., 106, 33951–33964,, 2001. a

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, 3 pp.,, 2017. a

Stevens, C. M., Verjans, V., Lundin, J. M., Kahle, E. C., Horlings, A. N., Horlings, B. I., and Waddington, E. D.: UWGlaciology/CommunityFirnModel: Version 1.0.5 of the Community Firn Model,, 2019. a

Sturm, M., Holmgren, J., König, M., and Morris, K.: The thermal conductivity of seasonal snow, J. Glaciol., 43, 26–41,, 1997. a

Tange, O.: GNU Parallel: the command-line power tool, login: The USENIX Magazine, 36, 42–47,, 2011. a

The IMBIE Team: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 558, 219–222,, 2018. a

Trudinger, C. M., Enting, I. G., Etheridge, D. M., Francey, R. J., Levchenko, V. A., Steele, L. P., Raynaud, D., and Arnaud, L.: Modeling air movement and bubble trapping in firn, J. Geophys. Res.-Atmos., 102, 6747–6763,, 1997. a

Tusima, K.: Grain Coarsening of Ice Particles Immersed in Pure Water, J. Jpn. Soc. Snow Ice, 40, 155–165,, 1978. a

Van Dusen, M. S.: Thermal conductivity of non-metallic solids, in: International critical tables of numerical data, physics, chemistry and technology, edited by: Washburn, E. W., 5, 216–217, McGraw-Hill, 1929. a

van Kampenhout, L., Lenaerts, J. T., 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, J. Adv. Model. Earth Sy., 9, 2583–2600,, 2017. a, b

van Wessem, J. M., van de Berg, W. J., Noël, B. P. Y., van Meijgaard, E., Amory, C., Birnbaum, G., Jakobs, C. L., Krüger, K., Lenaerts, J. T. M., Lhermitte, S., Ligtenberg, S. R. M., Medley, B., Reijmer, C. H., van Tricht, K., Trusel, L. D., van Ulft, L. H., Wouters, B., Wuite, J., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Antarctica (1979–2016), The Cryosphere, 12, 1479–1498,, 2018. a

Vandecrux, 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, J. Geophys. Res.-Ea. Surf. 123, 2563–2576,, 2018. a

Vandecrux, 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 compilation reveals widespread decrease of firn air content in western Greenland, The Cryosphere, 13, 845–859,, 2019. a

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., Simonsen, S. B., Heilig, A., Samimi, 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 Discuss.,, in review, 2020. a

Verjans, V., Leeson, A. A., Stevens, C. M., MacFerrin, M., Noël, B., and van den Broeke, M. R.: Development of physically based liquid water schemes for Greenland firn-densification models, The Cryosphere, 13, 1819–1842,, 2019. a, b, c, d, e, f

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, Geosci. Model Dev., 5, 773–791,, 2012. a, b

Voller, V. R., Swaminathan, C. R., and Thomas, B. G.: Fixed grid techniques for phase change problems: A review, Int. J. Numer. Meth. Eng., 30, 875–898,, 1990. a

Witrant, E., Martinerie, P., Hogan, C., Laube, J. C., Kawamura, K., Capron, E., Montzka, S. A., Dlugokencky, E. J., Etheridge, D., Blunier, T., and Sturges, W. T.: A new multi-gas constrained model of trace gas non-homogeneous transport in firn: evaluation and behaviour at eleven polar sites, Atmos. Chem. Phys., 12, 11465–11483,, 2012. a

Yen, Y.-C.: Review of thermal properties of snow, ice and sea ice, Tech. Rep. Report 81-10, Cold Regions Research and Engineering Lab, 1981. a

Short summary
Understanding processes in snow (firn), including compaction and airflow, is important for calculating how much mass the ice sheets are losing and for interpreting climate records from ice cores. We have developed the open-source Community Firn Model to simulate these processes. We used it to compare 13 different firn compaction equations and found that they do not agree within 10 %. We also show that including firn compaction in a firn-air model improves the match with data from ice cores.