The Community Firn Model (CFM) v1.0

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 Team, 2018). 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 Schwander, 2000). 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 Langway, 1980; Maeno and Ebinuma, 1983). The first zone extends from the surface, where density is often assumed to be ∼ 300–350, to 550 kg m⁻³. In zone 1, densification is usually considered to be due to grain-boundary sliding and settling (Alley, 1987). In zone 2, which spans the densities between 550 and ∼ 830 kg m⁻³, densification occurs due to sintering processes (Gow, 1975). Near a density of 830 kg m⁻³, 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) (Bader, 1954). 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 Whillans, 1994), 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 (Arzt, 1982) 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 Wingham, 2014), 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 Box, 2001; 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 1 List 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.

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 $\dot{b}$, and current density:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\text{d}\mathit{\rho }}{\text{d}t}}=f(T,\dot{b},\mathit{\rho }).\end{array}$$

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

$$\begin{array}{}\text{(2)}& {\mathit{\rho }}_{\text{new}}={\mathit{\rho }}_{\text{old}}+(\text{d}\mathit{\rho }/\text{d}t)\text{d}t.\end{array}$$

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 $\dot{b}$ by the relation $\mathit{\sigma }\left(z\right)=\dot{b}g\mathit{\tau }\left(z\right)$, 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 $\stackrel{\mathrm{‾}}{\dot{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 age_z,t is determined by the integrated accumulation-rate ($\dot{b}$) history (Li and Zwally, 2011, 2015):

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{\dot{b}}(z,t)={\displaystyle \frac{\mathrm{1}}{{\text{age}}_{z,t}}}\underset{t-{\text{age}}_{z,t}}{\overset{t}{\int }}\dot{b}\left(t^{\prime }\right)\text{d}{t}^{\prime }.\end{array}$$

The CFM uses $\stackrel{\mathrm{‾}}{\dot{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 Zwally, 2011), 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; Schytt, 1958; Robin, 1958), and (2) the firn's densification rate has an Arrhenius dependence on temperature. These assumptions can be incorporated into a densification-rate equation:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\text{d}\mathit{\rho }}{\text{d}t}}=c({\mathit{\rho }}_{\text{ice}}-\mathit{\rho }),\end{array}$$

with

$$\begin{array}{}\text{(5)}& c=k\exp \left(-{\displaystyle \frac{Q}{RT}}\right){\dot{b}}^a,\end{array}$$

where k and a are constants, Q is the Arrhenius activation energy (kJ mol⁻¹), R is the gas constant (8.314 $\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), and T is the temperature (K). For HL, c in Eq. (4) is given by

$$\begin{array}{}\text{(6)}& \begin{array}{rl}c& =c_{\mathrm{0}}=\mathrm{11}\exp \left(-{\displaystyle \frac{\mathrm{10.16}}{RT}}\right){\dot{b}}^{\mathrm{1.0}}(\mathit{\rho }\le \mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}})\\ & =c_{\mathrm{1}}=\mathrm{575}\exp \left(-{\displaystyle \frac{\mathrm{21.4}}{RT}}\right){\dot{b}}^{\mathrm{0.5}}(\mathit{\rho }>\mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}).\end{array}\end{array}$$

HL uses units of meters water equivalent per year (m w.e. a⁻¹) for $\dot{b}$. T in the original model was the mean annual site temperature T_m, 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 $\dot{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⁻³), which are numerically equivalent to units grams per cubic centimeter (g cm⁻³).

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 ρ>550 kg m⁻³, the densification rate is given by

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\text{d}\mathit{\rho }}{\text{d}t}}={\mathit{\rho }}_{\text{i}}{A}_{\mathrm{0}}\text{exp}\left({\displaystyle \frac{-Q}{RT}}\right)f{\mathit{\sigma }}_{\text{d}}ef{f}^n,\end{array}$$

where $A_{\mathrm{0}}=\mathrm{2.54}\times {\mathrm{10}}^{\mathrm{4}}$ ${\mathrm{MPa}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$, the activation energy Q is 60 kJ mol⁻¹, σ_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

$$\begin{array}{}\text{(8)}& f={\mathrm{10}}^{\mathit{\alpha }{\mathit{\rho }}^{\mathrm{3}}+\mathit{\beta }{\mathit{\rho }}^{\mathrm{2}}+\mathit{\delta }\mathit{\rho }+\mathit{\gamma }}(\mathrm{550}\le \mathit{\rho }\le \mathrm{800}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}),\end{array}$$

with $\mathit{\alpha }=-\mathrm{37.455}$, β=99.743, $\mathit{\delta }=-\mathrm{95.027}$, and γ=30.673.

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

$$\begin{array}{}\text{(9)}& f={\displaystyle \frac{\mathrm{3}}{\mathrm{16}}}(\mathrm{1}-\mathit{\rho }/{\mathit{\rho }}_i)/(\mathrm{1}-(\mathrm{1}-\mathit{\rho }/{\mathit{\rho }}_i{)}^{\mathrm{1}/\mathrm{3}}{)}^{\mathrm{3}}(\mathit{\rho }>\mathrm{800}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}).\end{array}$$

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

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\text{d}D}{\text{d}t}}=\mathit{\gamma }\left({\displaystyle \frac{P}{D^{\mathrm{2}}}}\right)\left(\mathrm{1}-{\displaystyle \frac{\mathrm{5}}{\mathrm{3}}}D\right)(\mathit{\rho }\le \mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}),\end{array}$$

where $D=\mathit{\rho }/{\mathit{\rho }}_{\text{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):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& {\displaystyle \frac{\text{d}D}{\text{d}t}}=\\ & \mathrm{4.1817}\times {\mathrm{10}}^{\mathrm{4}}\exp \left(-{\displaystyle \frac{E_{\text{A}}}{RT}}\right){\left(D^{\mathrm{2}}{D}_{\mathrm{0}}\right)}^{\mathrm{1}/\mathrm{3}}{\left({\displaystyle \frac{a}{\mathit{\pi }}}\right)}^{\mathrm{1}/\mathrm{2}}{\left({\displaystyle \frac{\mathrm{4}\mathit{\pi }P}{\mathrm{3}aZD}}\right)}^{\mathrm{3}}\\ & (\mathit{\rho }>\mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}),\end{array}\end{array}$$

where E_A is the activation energy given as 60 kJ mol⁻¹, 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. D₀ 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 T_m (K), given by

$$\begin{array}{}\text{(12)}& D_{\mathrm{0}}=\mathrm{0.00226}{T}_{\text{m}}+\mathrm{0.03}.\end{array}$$

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

Buizert et al. (2015) described an issue in implementing the Goujon et al. (2003) model, which we review here. In the event that D₀≥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=D₀, 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 D₀ given by Eq. (12) to a maximum value of 0.59, which occurs for temperatures greater than $\sim -\mathrm{25}$ ^∘C. This value of D₀ corresponds to a density of 541 kg m⁻³, 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⁻³. We follow the suggestion in Buizert et al. (2015) and put the zone 1–zone 2 transition at relative density $D_{\mathrm{0}}^{\prime }=D_{\mathrm{0}}+\mathit{ϵ}$, where ϵ is a small number. The densification rate in zone 2 is still calculated using Eq. (11) using D₀ 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 Zwally, 2002, 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:

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\text{d}\mathit{\rho }}{\text{d}t}}=\mathit{\beta }\mathrm{8.36}(\mathrm{273.2}-T_{\text{K}}{)}^{-\mathrm{2.061}}\stackrel{\mathrm{‾}}{\dot{b}}({\mathit{\rho }}_{\text{i}}-\mathit{\rho }),\end{array}$$

where T_K is the firn temperature as a function of time and depth in Kelvin, and $\stackrel{\mathrm{‾}}{\dot{b}}$ is the mean accumulation rate over the lifetime of a parcel of firn (Eq. 3) (m w.e. a⁻¹). The difference between LZ11 and LZ15 is in the parameter β. For LZ11,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}\mathit{\beta }=& {\mathit{\beta }}_{\mathrm{1}}=-\mathrm{9.788}+\mathrm{8.996}{\dot{b}}_{\text{m}}-\mathrm{0.6165}{T}_{\text{m,c}}\\ & (\mathit{\rho }\le \mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}})\\ =& {\mathit{\beta }}_{\mathrm{2}}={\mathit{\beta }}_{\mathrm{1}}/(-\mathrm{2.0178}+\mathrm{8.4043}{\dot{b}}_{\text{m}}-\mathrm{0.0932}{T}_{\text{m,c}})\\ & (\mathit{\rho }>\mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}).\end{array}\end{array}$$

For LZ15,

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\mathit{\beta }=& {\mathit{\beta }}_{\mathrm{1}}=-\mathrm{1.218}-\mathrm{0.403}{T}_{\text{m,c}}\\ & (\mathit{\rho }\le \mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}})\\ =& {\mathit{\beta }}_{\mathrm{2}}={\mathit{\beta }}_{\mathrm{1}}(\mathrm{0.792}-\mathrm{1.080}{\dot{b}}_{\text{m}}+\mathrm{0.00465}{T}_{\text{m,c}})\\ & (\mathit{\rho }>\mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}),\end{array}\end{array}$$

where T_m,c is the mean annual surface temperature in Celsius and ${\dot{b}}_{\text{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

$$\begin{array}{}\text{(16)}& \mathit{\beta }={\mathit{\beta }}_{\mathrm{1}}={\mathit{\beta }}_{\mathrm{2}}=\mathrm{76.138}-\mathrm{0.28965}{T}_{\text{m}}.\end{array}$$

The mean annual surface temperature T_m 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):

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& {\displaystyle \frac{\text{d}\mathit{\rho }}{\text{d}t}}=k_{\text{c}}({\mathit{\rho }}_{\text{i}}-\mathit{\rho })\exp (-E_{\text{c}}/RT)\mathit{\sigma }/r^{\mathrm{2}}\\ & {\displaystyle \frac{\text{d}{r}^{\mathrm{2}}}{\text{d}t}}=k_{\text{g}}\exp (-E_{\text{g}}/RT)\\ & {\displaystyle \frac{\text{d}\mathit{\sigma }}{\text{d}t}}=\dot{b}g,\end{array}\end{array}$$

where k_c and k_g are empirically derived constants. The water-molecule diffusion activation energy E_c is 60 kJ mol⁻¹, and the grain-growth activation energy E_g is 42.4 kJ mol⁻¹. 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

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& c_{\mathrm{0}}=\mathrm{0.07}\dot{b}g\exp \left(-{\displaystyle \frac{E_{\text{c}}}{RT}}+{\displaystyle \frac{E_{\text{g}}}{R{T}_{\text{m}}}}\right)(\mathit{\rho }\le \mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}})\\ & c_{\mathrm{1}}=\mathrm{0.03}\dot{b}g\exp \left(-{\displaystyle \frac{E_{\text{c}}}{RT}}+{\displaystyle \frac{E_{\text{g}}}{R{T}_{\text{m}}}}\right),(\mathit{\rho }>\mathrm{550}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}),\end{array}\end{array}$$

with $\dot{b}$ the mass accumulation rate ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$) and T_m 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 c₀ and c₁ in Eq. (18) are multiplied by additional tuning coefficients. For LIG,

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{c}_{\mathrm{0}}^{\text{(LIG)}}& =[\mathrm{1.435}-\mathrm{0.151}\ln (\dot{b}\left)\right]c_{\mathrm{0}}^{\text{(ART-S)}}\\ c_{\mathrm{1}}^{\text{(LIG)}}& =[\mathrm{2.366}-\mathrm{0.293}\ln (\dot{b}\left)\right]c_{\mathrm{1}}^{\text{(ART-S)}},\end{array}\end{array}$$

and for KM,

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{c}_{\mathrm{0}}^{\text{(KM)}}& =[\mathrm{1.042}-\mathrm{0.0916}\ln (\dot{b}\left)\right]c_{\mathrm{0}}^{\text{(ART-S)}}\\ c_{\mathrm{1}}^{\text{(KM)}}& =[\mathrm{1.734}-\mathrm{0.2039}\ln (\dot{b}\left)\right]c_{\mathrm{1}}^{\text{(ART-S)}}.\end{array}\end{array}$$

For LIG and KM, $\dot{b}$ has units of kilograms per square meter per year ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{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 T_m in the grain-growth term of the Arrhenius factor; Steger et al. (2017) modified this to use T_m, 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 c₀ and c₁ in Eq. (18) by additional factors:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{c}_{\mathrm{0}}^{\text{(SIM)}}& =f_{\mathrm{0}}{c}_{\mathrm{0}}^{\text{(ART-S)}}\\ c_{\mathrm{1}}^{\text{(SIM)}}& =f_{\mathrm{1}}{\displaystyle \frac{\mathrm{61.7}}{{\dot{b}}^{\mathrm{0.5}}}}\exp \left({\displaystyle \frac{-\mathrm{3800}}{R{T}_m}}\right)c_{\mathrm{1}}^{\text{(ART-S)}}.\end{array}\end{array}$$

Simonsen et al. (2013) found that the parameters f₀ and f₁ 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 f₀=0.8 and f₁=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):

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& {\displaystyle \frac{\text{d}\mathit{\rho }}{\text{d}t}}={\displaystyle \frac{\mathit{\rho }\mathit{\sigma }}{\mathit{\eta }}}\\ & \mathit{\eta }=f_{\mathrm{1}}{f}_{\mathrm{2}}{\mathit{\eta }}_{\mathrm{0}}{\displaystyle \frac{\mathit{\rho }}{c_{\mathit{\eta }}}}\exp \left(a_{\mathit{\eta }}(T_{\text{melt}}-T)+b_{\mathit{\eta }}\mathit{\rho }\right),\end{array}\end{array}$$

where ρ is the density, T_melt is the melting temperature, ${\mathit{\eta }}_{\mathrm{0}}=\mathrm{7.62237}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$, $a_{\mathit{\eta }}=\mathrm{0.1}{\mathrm{K}}^{-\mathrm{1}}$, $b_{\mathit{\eta }}=\mathrm{0.023}{\mathrm{m}}^{\mathrm{3}}{\mathrm{kg}}^{-\mathrm{1}}$. f₁ and f₂ are snow-viscosity correction factors. f₁ 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. f₂ accounts for angular grains. Following Langen et al. (2017) and van Kampenhout et al. (2017), f₂ is set to 4 in the CFM. Vionnet and others (2012) give $c_{\mathit{\eta }}=\mathrm{250}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, but also following van Kampenhout et al. (2017), we set $c_{\mathit{\eta }}=\mathrm{358}\mathrm{kg}{\mathrm{m}}^{-\mathrm{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⁻³:

$$\begin{array}{}\text{(23)}& \dot{\mathit{ϵ}}={\displaystyle \frac{k_{\mathrm{0}}^{*}}{{\mathit{\rho }}_{\text{w}}g}}\left({\displaystyle \frac{{\mathit{\rho }}_{\text{i}}-\mathit{\rho }}{\mathit{\rho }}}\right)\left(\mathrm{1}-{\stackrel{\mathrm{‾}}{M}}_{\mathrm{0}}m\right){\displaystyle \frac{\mathrm{1}}{H\left(\mathit{\tau }\right)}}\exp \left(-{\displaystyle \frac{E_{\mathit{\alpha }}}{RT}}\right)\mathit{\sigma },\end{array}$$

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

$$\begin{array}{}\text{(24)}& H\left(\mathit{\tau }\right)=\underset{{\mathit{\tau }}_{\mathrm{0}}}{\overset{\mathit{\tau }}{\int }}\exp \left(-{\displaystyle \frac{E_{\mathit{\alpha }}}{RT\left({\mathit{\tau }}^{\prime }\right)}}\right)\text{d}{\mathit{\tau }}^{\prime }.\end{array}$$

In the above equations, ρ_i and ρ_w are the densities of ice and water, respectively, $k_{\mathrm{0}}^{*}$ is a densification constant, and τ is the age of a parcel of firn. τ₀ 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 ${\stackrel{\mathrm{‾}}{M}}_{\mathrm{0}}$ is a scaling constant. The CFM by default uses the preferred activation energy E_α presented in Morris and Wingham (2014) of 110 kg m⁻³, though the authors also test their model with values of 60 and 200 kJ mol⁻¹. 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 $k_{\mathrm{0}}^{*}=\mathrm{11}$ $\mathrm{m}\mathrm{w}.\mathrm{e}{.}^{-\mathrm{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_α−E_H (revising Table 2 from Morris and Wingham, 2014), which can be used to calculate the corrected value of $k_{\mathrm{0}}^{*}$ with the following equation:

$$\begin{array}{}\text{(25)}& k_{\mathrm{0}}^{*}=\mathrm{11}\exp \left(-{\displaystyle \frac{E^{*}-(E_{\mathit{\alpha }}-E_H)}{R{T}_{\text{m}}}}\right),\end{array}$$

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

Table 2 E_α−E_H (kJ mol⁻¹), revised from Table 2 of Morris and Wingham (2014).

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 (Patankar, 1980). 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:

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\text{D}\mathit{\delta }}{\text{D}t}}=\mathrm{\Omega }(t,z){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\delta }}{\partial z^{\mathrm{2}}}},\end{array}$$

where $\frac{\text{D}\mathit{\delta }}{\text{D}t}$ 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 δ¹⁸O 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):

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{\displaystyle \frac{\partial C}{\partial t}}& ={\displaystyle \frac{\mathrm{1}}{{\mathit{\varphi }}_{\text{op}}}}{\displaystyle \frac{\partial }{\partial z}}\left[{\mathit{\varphi }}_{\text{op}}(t,z){\mathit{\kappa }}_{\text{eff}}(t,z)\left({\displaystyle \frac{\partial C}{\partial z}}-{\displaystyle \frac{\mathrm{\Delta }mg}{RT}}+\mathrm{\Omega }{\displaystyle \frac{\partial T}{\partial z}}\right)\right]\\ & -w_{\text{air}}\left(z\right){\displaystyle \frac{\partial C}{\partial z}},\end{array}\end{array}$$

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⁻¹); Ω is the thermal diffusion sensitivity (K⁻¹), which is specific to individual gases (Severinghaus et al., 2001); w_air(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): ${\mathit{\kappa }}_{\text{eff}}={\mathit{\kappa }}_{\text{FA}}/\mathit{\tau }$. 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 w_air(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 ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\text{co}}$, the CFM uses the relationship published by Martinerie et al. (1994):

$$\begin{array}{}\text{(28)}& {\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\text{co}}={\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\text{ice}}}}+\mathrm{6.95}\times {\mathrm{10}}^{-\mathrm{7}}{T}_{\text{m}}-\mathrm{4.3}\times {\mathrm{10}}^{-\mathrm{5}}\right)}^{-\mathrm{1}},\end{array}$$

with temperature in Kelvin (K) and ρ in kilograms per cubic meter (kg m⁻³). 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⁻³. 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 ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\text{co}}$ and its depth z_co slightly deeper than $z\left({\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\text{co}}\right)$.

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):

$$\begin{array}{}\text{(29)}& {\mathit{\varphi }}_{\text{cl}}=\mathrm{0.37}{\mathit{\varphi }}_{\text{tot}}{\left({\displaystyle \frac{{\mathit{\varphi }}_{\text{tot}}}{{\mathit{\varphi }}_{\text{co}}}}\right)}^{-\mathrm{7.6}}.\end{array}$$

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⁻³ from ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\text{co}}$ (Blunier and Schwander, 2000).

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 T_melt), 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 Lesaffre, 1998) 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⁻³ (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):

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{\text{d}r}{\text{d}t}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}r}}{k}_{\text{g}}\exp (-E_{\text{g}}/RT),\end{array}$$

with grain-growth activation energy $E_{\text{g}}=\mathrm{42.4}\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$ and constant $k_{\text{g}}=\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹ 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

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{\displaystyle \frac{\text{d}r}{\text{d}t}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{8}\times {\mathrm{10}}^{\mathrm{9}}{r}^{\mathrm{2}}}}\times min\left[{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}\left(\mathrm{1.28}\times {\mathrm{10}}^{-\mathrm{8}}+\mathrm{4.22}\right.\right.\\ & \left.\left.\times {\mathrm{10}}^{-\mathrm{10}}{\mathit{\theta }}^{\mathrm{3}}\right),\mathrm{6.94}\times {\mathrm{10}}^{-\mathrm{8}}\right],\end{array}\end{array}$$

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 r₀=0.1 mm. For Eq. (31), the CFM uses the empirical formula for surface grain size as a function of T_m (^∘C) and $\dot{b}$ (m w.e. a⁻¹) given by Linow et al. (2012):

$$\begin{array}{}\text{(32)}& r_{\mathrm{0}}=b_{\mathrm{0}}+b_{\mathrm{1}}{T}_{\text{m}}+b_{\mathrm{2}}\dot{b},\end{array}$$

with constants b₀=0.781, b₁=0.0085, and $b_{\mathrm{2}}=-\mathrm{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.4 ^∘C, and the annual accumulation averages 0.23 m ice eq. a⁻¹. 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⁻³. The surface density was held constant at 300 kg m⁻³, 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×1 m² firn column (m³ m⁻²) above depth z, given by

$$\begin{array}{}\text{(33)}& \text{DIP}\left(z\right)=\underset{\mathrm{0}}{\overset{z}{\int }}\mathit{\varphi }\left(z^{\prime }\right)\text{d}{z}^{\prime }=\underset{\mathrm{0}}{\overset{z}{\int }}{\displaystyle \frac{{\mathit{\rho }}_{\text{i}}-\mathit{\rho }\left(z^{\prime }\right)}{{\mathit{\rho }}_{\text{i}}}}\text{d}{z}^{\prime },\end{array}$$

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⁻³. 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 DIP₁₅, DIP₈₀, and DIP_tot, 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 $\dot{b}$, firn compaction rate v_fc, vertical ice velocity v_ice at the bottom of the firn column due to dynamic ice-sheet processes, and vertical bedrock motion rate v_bed:

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\text{d}H}{\text{d}t}}=\dot{b}+v_{\text{fc}}+v_{\text{ice}}+v_{\text{bed}}.\end{array}$$

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 v_ice 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 v_bed 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⁻³, but in reality, the BCO depth does not correspond to an exact density. Nevertheless, here we use the 830 kg m⁻³ 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 1 Profiles 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 DIP₁₅, DIP₈₀, and DIP_tot; the depth and age of the 830 kg m⁻³ 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' DIP_tot is not a direct comparison metric (as opposed to comparing DIP₁₅). 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⁻³ at 126 m and 917 kg m⁻³ at 180 m; the DIP in this interval is 0.07 m, or 0.3 % of DIP_tot. We thus consider DIP_tot to be a worthy metric for comparison.

Near the surface (zone 1 densification), the models show a variety of responses. The mean DIP₁₅ is 7.75 m with σ=0.81 m, 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 DIP₁₅ 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=\mathrm{0}\mathrm{m})=\mathrm{0}$ m and ρ_surface=300 kg m⁻³). 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 DIP₈₀, and ART-S predicts the lowest DIP₈₀. The mean of the models' DIP₈₀ 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 DIP₈₀ is similar to the other models, but its density is ∼100 kg m⁻³ less than the other models. CRO's DIP_tot 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⁻³ density horizon (excluding CRO) is 67 m, and σ=9.07 m (CV =13 %). The mean age is 224 years with σ=29 years (CV =13 %). In general, models that predict deeper DEP₈₃₀ also predict older AGE₈₃₀ when compared with the others, but this is not universally true. For example, LZ15 predicts DEP₈₃₀=65.46 m and AGE₈₃₀=232 years. LIG predicts a similar DEP₈₃₀ of 64.37 m but a much younger AGE₈₃₀ of 214 years (1 m of firn at Summit at this density is roughly 4 years of accumulation). Likewise, BAR predicts DEP₈₃₀=70.03 m and AGE₈₃₀=231 years, while HEL predicts a deeper DEP₈₃₀ (73.09 m) and younger AGE₈₃₀ (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 $\dot{b}$ sensitivity. In this case, these faulty sensitivities compensate for each other; the model will predict dρ∕dt and AGE₈₃₀ reasonably well. However, due to the overestimated $\dot{b}$ sensitivity, not enough material has been accumulated over the time required to reach a density of 830 kg m⁻³, and the model will predict a DEP₈₃₀ that is too shallow.

Table 3 Model 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⁻³ 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.

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⁻¹) for each model from October in 2008–2016, 2016–2018, and 2008–2016 (columns dH∕dt_08–16, dH∕dt_16–18, and dH∕dt_08–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.042 cm a⁻¹, 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.13 cm a⁻¹) and MW the smallest (2.61 cm a⁻¹). The mean of the models' trends for these two years is 3.22 cm a⁻¹ with σ=0.48 cm a⁻¹ (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⁻¹. 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⁻¹, and σ=0.27 cm a⁻¹ (CV =32 %). ART-T predicts the smallest trend in dH∕dt (0.28 cm a⁻¹; 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.19 cm a⁻¹), 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 DIP₁₅, 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 2 Model-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 (Alley, 2004; Cuffey and Clow, 1997) 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⁻¹ and mean annual temperature −47.5 ^∘C.

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 Brook, 1999) 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 δ¹⁵N and δ⁴⁰Ar isotopes; δ¹⁵N and δ⁴⁰Ar∕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 δ¹⁵N and δ⁴⁰Ar in the firn at Summit. We forced the CFM with temperature and accumulation-rate histories from the GISP2 ice core (Cuffey and Clow, 1997; Alley, 2000, 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⁻¹ and a temperature of −47.5 ^∘C, 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 4 Measured and modeled δ¹⁵N and δ⁴⁰Ar∕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 δ¹⁵N 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 δ¹⁵N and δ⁴⁰Ar∕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 δ¹⁵N 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 5 As 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 δ¹⁵N and δ⁴⁰Ar∕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.5 ^∘C and 0.07 m ice eq. a⁻¹, 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⁻¹ 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 δ¹⁵N and δ⁴⁰Ar∕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 δ¹⁵N values. The transient model matches the observed magnitude and rate of the δ¹⁵N and δ⁴⁰Ar∕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 Langway, 1980) 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.